June 11, 2008 

Preprint typeset using IAT^X style cmulatcapj v. 03/07/07 



00 
O 
O 
(N 

C 
3 



43 
Oh- 

6 



o 

O 

oo 
O 



X 



EVOLUTION OF MIGRATING PLANETS UNDERGOING GAS ACCRETION^ 

Gennaro D'Angelo 1 

NASA Ames Research Center, Space Science and Astrobiology Division, MS 245-3, Moffett Field, CA 94035 

AND 

Stephen H. Lubow 

Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218 

June 11, 2008 

ABSTRACT 

We analyze the orbital and mass evolution of planets that undergo run-away gas accretion by 
means of two- and three-dimensional hydrodynamic simulations. The disk torque distribution per 
unit disk mass as a function of radius provides an important diagnostic for the nature of the disk- 
planet interactions. We first consider torque distributions for nonmigrating planets of fixed mass and 
show that there is general agreement with the expectations of resonance theory. We then present 
results of simulations for mass-gaining, migrating planets. For planets with an initial mass of 5 Earth 
masses (Me), which are embedded in disks with standard parameters and which undergo run-away gas 
accretion to one Jupiter mass (Mj), the torque distributions per unit disk mass are largely unaffected 
by migration and accretion for a given planet mass. The migration rates for these planets are in 
agreement with the predictions of the standard theory for planet migration (Type I and Type II 
migration). The planet mass growth occurs through gas capture within the planet's Bondi radius at 
lower planet masses, the Hill radius at intermediate planet masses, and through reduced accretion at 
higher planet masses due to gap formation. During run-away mass growth, a planet migrates inwards 
by only about 20% in radius before achieving a mass of ~ lMj. For the above models, we find no 
evidence of fast migration driven by coorbital torques, known as Type III migration. We do find 
evidence of Type III migration for a fixed mass planet of Saturn's mass that is immersed in a cold and 
massive disk. In this case the planet migration is assumed to begin before gap formation completes. 
The migration is understood through a model in which the torque is due to an asymmetry in density 
between trapped gas on the leading side of the planet and ambient gas on the trailing side of the 
planet. 

Subject headings: accretion, accretion disks — hydrodynamics — methods: numerical — planetary 
systems: formation — planetary systems: protoplanctary disks — solar system: 
formation 



1. INTRODUCTION 

In the core accretion picture of planet formation (Bo- 
denheimer & Pollack 1986; Wuchterl 1991; Pollack et al. 
1996; Hubickyj et al. 2005, and references therein), a 
small mass solid core initially rapidly accretes solid ma- 
terial, followed by a slow evolution phase of gas and 
solid accretion. During this slow evolution phase, the 
planet is limited in its ability to accrete gas by the ther- 
mal heating caused by the impacting solids. Once the 
planet's gas mass is greater than its solid mass, typi- 
cally at several Earth masses, the planet undergoes "run- 
away" gas accretion, in which it can accrete whatever 
mass is provided to it. These processes have been treated 
by one-dimensional, spherically symmetric structure cal- 
culations in the above papers. 

On the other hand, multi-dimensional hydrodynami- 
cal calculations of a protostcllar disk interacting with 
the planet has revealed various flow properties of the 
gas, including the gap opening by tidal effects, previ- 
ously anticipated by one-dimensional disk models (Lin 
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& Papaloizou 1986). In addition, planet migration that 
results from disk-planet interactions has been analyzed 
by means of such simulations. Good agreement is often, 
but not always, found between the simulations and the 
expectations of theory (Nelson et al. 2000; Bate et al. 
2003; D'Angelo et al. 2003; Nelson & Benz 2003; Li et al. 
2005; D'Angelo et al. 2006). These calculations typically 
do not include the mass evolution of the planet. Usu- 
ally they apply accretion boundary conditions onto the 
planet as a means of modelling the run-away gas accre- 
tion process. One aim of this paper is to analyze the 
effects of planet mass growth on migration. 

Several controversies remain on the effects of gas. The 
role of coorbital torques on planet migration, in the sub- 
giant mass range, is not well understood. Masset & Pa- 
paloizou (2003, hereafter MP03) suggested on the basis 
of a model and simulations that a fast mode of migration 
(sometimes called Type III migration) can occur due to 
strong coorbital torques. Ogilvie & Lubow (2006, here- 
after OL06) found support for the concept of coorbital 
dominated migration under certain conditions. At higher 
grid resolution under the conditions specified by MP03, 
simulations by D'Angelo et al. (2005, hereafter DBL05) 
found that the migration rate was much slower. 

Another subject of interest is how planet masses may 
be limited by a reduction in the gas accretion rate. Lin 
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& Papaloizou (1986) proposed such a reduction by tidal 
torques that open a gap about the orbit of the planet. 
The value of the highest planet mass achieved in the 
presence of gap opening is somewhat controversial. Some 
studies (Lubow et al. 1999; Bate et al. 2003; D'Angelo 
et al. 2003) have suggested that the maximum planet 
mass is about 6-10 Mj, corresponding to the upper limit 
of the observed range of extrasolar planets (Marcy et al. 
2005; Butler et al. 2006). This limit suggests that some 
other process, such as disk dispersal or other self-limiting 
feedback on planetary accretion, is responsible for the 
lower masses (~ 1 Mj) typically found observationally. 
Other studies suggest that the tidal limit is ~ 1 Mj and 
therefore no additional process is required to explain the 
typical masses (e.g., Dobbs-Dixon et al. 2007). 

We will address these and other issues in this paper by 
analyzing the orbital evolution of a mass-gaining planet 
embedded in a gas disk. In section 2 we analyze the 
torque distributions for planets of constant mass on fixed 
circular orbits. In section 3 we analyze the orbital and 
mass evolution of migrating planets that undergo run- 
away mass accretion. Section 4 describes a model that 
appears to exhibit migration that is dominated by coor- 
bital torques, i.e., Type III migration. Section 5 contains 
the summary and discussion. 

2. TORQUE DISTRIBUTION FOR A 
NON-MIGRATING PLANET 

Disk-planet gravitational torques result in planet mi- 
gration (Goldreich & Tremaine 1980; Lin & Papaloizou 
1993; Ward 1997). The distribution of torque with disk 
radius provides a means of connecting the theory with 
simulations. In this section, we model the disk as a three- 
dimensional system and consider fixed mass planets on 
fixed circular orbits. The torque per unit radius for a 
planet embedded in a disk was previously considered in 
Bate et al. (2003). Here we reconsider the analysis with 
higher resolution, especially in the coorbital region, and 
apply the torque distribution per unit disk mass. 

2.1. Numerical Procedure 

In this section, we describe the torques exerted by a 
disk on an embedded planet with mass, M p , equal to 
1M E , 10 M e , 0.3 Mj, and lMj. For the two small- 
est mass planets we consider, the planet's Hill radius is 
smaller than the vertical disk thickness of several per- 
cent of the distance to the star. For the two largest mass 
planets, the Hill radius is comparable or larger than the 
disk thickness. 

2.1.1. Disk Model 

We use spherical polar coordinates {R, 8, <fi}, with the 
origin located at the star-planet center of mass. The 
reference frame corotates with the star-planet system. 
The planet's orbit lies in the plane 8 = tt/2. The disk 
is assumed to be symmetric with respect to this plane, 
hence only the disk's northern hemisphere (i.e., the vol- 
ume 8 < tt/2) is simulated. 

We assume that the material in the disk is locally 
isothermal and that the pressure p is given by 

p{R,e,<j>) = p{R,0A)c 2 s {r), (1) 
where p(R, 8, <j)) is the mass density. Quantity c s (r) is the 
gas sound speed, which is taken to be a function of cylin- 
drical radius r = Rsmd. The aspect ratio of the disk, 



H/r, is taken to be constant and equal to 0.05. There- 
fore, the temperature distribution in the disk is only a 
function of the distance from the disk's rotation axis, r, 
and decreases as c 2 oc 1/r. Viscous forces are calculated 
by adopting the stress tensor for a Newtonian fluid (Mi- 
halas & Weibel Mihalas 1999) with constant kinematic 
viscosity, v and zero bulk viscosity. Disk self-gravity is 
ignored. In Appendix C, we discuss some effects of disk 
self-gravity and of the axisymmetric component of disk 
gravity on the migration rates. 

2.1.2. Disk and Planet Parameters 

We adopt the stellar mass M s as unit of mass, 
the orbital radius a as unit of length, and flp 1 = 

1/2 

[G (M s + M p )/a 3 ] as unit of time. In converting 
to dimensional units we consider a = 5.2 AU and M s = 
1M . 

The disk extends from to 2tt in azimuth around the 
star and, in radius, from 0.4 to cither 4.0 (Jupiter-mass 
case) or 2.5 (lower mass cases). In the ^-direction, the 
disk domain extends above the midplane (8 = tt/2) for 
10 degrees, comprising 3.5 pressure scale heights, H. The 
initial mass density distribution is independent of <f> 1 has 
a Gaussian profile in the ^-direction, and has a radial 
profile proportional to i?~ 3 / 2 , so that the initial (un- 
perturbed) surface density varies as R- 1 / 2 . We adopt 
a constant dimensionless kinematic viscosity v equal to 
10 -5 , corresponding to a turbulent viscosity parameter 
a = 0.004 at the cylindrical radius r = 1 (5.2 AU). 

As mentioned above, we perform calculations for four 
planet masses: M p = 3 x 10~ 6 , 3 x 10~ 5 , 3 x 10 -4 , and 
1 x 10 -3 , which correspond, respectively, to 1 Me, 10 Me, 
0.3 Mj, and 1 Mj. The gravitational potential, $ p , of the 
planet is smoothed over a length e equal to 0.1 i?n and 
is given by 



where S is the distance from the planet and i?u is the 
Hill radius of the planet. 

2.1.3. Numerical Method 

The mass and momentum equations that describe the 
evolution of the disk (e.g., DBL05) are solved numer- 
ically by means of a finite-difference scheme that ap- 
plies an operator splitting procedure to perform the spa- 
tial integration of advection and source terms (Ziegler 
& Yorke 1997). The algorithm is second-order accurate 
in space and semi-sccond-order in time. The equations 
are discrctized over a mesh with constant grid spacing 
in each coordinate direction. Nested grids are used to 
enhance the numerical resolution in (arbitrarily large) 
regions around the planet (D'Angelo et al. 2002, 2003). 
This strategy allows the volume resolution to be in- 
creased by a factor 2 3 for each added grid level. These 
calculations are executed with grid systems involving 
5 levels of grid nesting. The linear base resolution is 
AR = a Ad = aA4> = 0.014 a. The linear resolution 
achieved in the coorbital region around the planet is ap- 
proximately 9 x 10 -4 a, which corresponds to ~ 0.01 i?H 
and ~ 0.1 i?H in the Jupiter-mass and Earth-mass cases, 
respectively. To quantify resolution effects in the Earth- 
mass case, we also applied a linear resolution twice as 
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high throughout the entire grid system (base resolution 
of 7x 10~ 3 a and resolution in the coorbital region around 
the planet of 4 x 10~ 4 a). The torques at the two reso- 
lutions, integrated over the disk domain, differ by about 
5%. 

The boundary condition near the planet involves re- 
moving gas from ~ 0.1 i?H of the planet at each timestep. 
The procedure for mass removal is described in more de- 
tail in section 3.1.1. In the calculations reported in sec- 
tion 2.3, the removed mass is not added to the planet's 
mass in order to keep it fixed. In sections 3 and 4 (as 
well as in Appendix A and C), we will present cases in 
which the planet's mass is augmented by the mass of the 
gas removed from the disk. 

The outer boundary of the disk domain is closed to 
both inflow and outflow, whereas the inner boundary al- 
lows outflow (material can flow out of the grid domain) 
but not inflow. Reflective and symmetry boundary con- 
ditions are applied at colatitude 6 = 6* m i n and at the disk 
mid-plane (6 = tt/2), respectively. 

Simulations are run for about 100 orbital periods. In 
models with 0.3 Mj and 1 Mj mass planets, the initial 
density distribution includes a gap along the planet's or- 
bit to account for an approximate balance between vis- 
cous and tidal torques, which reduces the relaxation time 
towards steady state. In all calculations discussed here, 
the flow achieves a fairly steady state within ~ 100 or- 
bits. 

2.2. Theoretical Considerations 
2.2.1. Torque Density 

Consider a cylindrical coordinate system {r, (j>, z} cen- 
tered on the star-planet center of mass. The disk torque 
along the rotation axis per unit radius exerted on the 
planet is given by 



alT 
dr 



(r,t) 



dzp(r,t)fy* p (r,t)), (3) 



where (X(t)) denotes the time-average of X over an orbit 
period centered about time t, p is the gas density, and 
$ p is the potential due to the planet (eq. 2). 

2.2.2. Radial Overlap Regions 

The linear theory of Lindblad resonances for disk- 
planet interactions demonstrates that the strongest con- 
tributing resonances have azimuthal wavenumbers m ~ 
r/H. This estimate comes from considering the so-called 
torque cutoff effect that arises from Lindblad resonances 
that lie close to the planet (Goldrcich & Trcmaine 1980; 
Ward 1986; Artymowicz 1993). As a consequence of the 
resonance condition, we expect the peak torque density 
to be at a distance of roughly H from the planet. The 
torque cutoff is not sharp and there are torque contribu- 
tions from resonances that lie closer than distance ~ H 
from the planet, although at a decreasing level as they 
get closer to the planet. As we will see, the numerical 
results show the torque density peak to be close to dis- 
tance H from the planet. However, the torque cutoff 
calculations assume that the orbits are such that the gas 
azimuthally passes by the planet, i.e, lies on circulating 
orbits. On the other hand, close to the planet's orbit, 
this assumption breaks down and the gas flows on librat- 
ing streamlines of the horseshoe orbit region. This region 



generally extends in the radial direction to a distance of 
about 3i?H from the planet's orbital radius, where i?H 
is the planet's Hill radius. But, close to the planet, the 
region becomes less extended radially, spanning only to 
approximately i?n- That is, the noncoorbital (circulat- 
ing) streamlines pass closest to the planet at a distance 
about equal to i?H (see streamline a in Lubow et al. 1999 
and Figure 5 in Bate et al. 2003). In the horseshoe orbit 
region, the corotational resonance can play a role. 

These two regions, the coorbital region (extending up 
to about 3 i?H from the planet's orbital radius) and Lind- 
blad torque region (extending beyond about distance 
H from the planet's orbital radius), overlap in a one- 
dimensional radial sense for planet-to-star mass ratios 
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1 fH 



(4) 



This condition does not necessarily imply a physical over- 
lap in two or three dimensions. But it does affect our 
interpretation of the torque density reduced to one di- 
mension, dT(r)/dr. The reason is that for a given radius 
r such that i?n < \r — o| < 3 i?H, the gas lies in either the 
coorbital (librating) or noncoorbital (circulating) region, 
depending on the azimuth. 

For the disk parameters considered in this section, the 
one-dimensional overlap occurs for planet masses greater 
than about 4.6 Me, which covers all, but one, of the 
planet masses considered. For a 1 Mj planet, this overlap 
occurs out to a radius of about 1.2 a or a radial distance 
of about AH from planet. 

The two regions physically overlap in a two- or three- 
dimensional sense, when the closest approach of all non- 
coorbital (circulating) streamlines, which occurs at a dis- 
tance ~ i?H from the planet, is greater than the distance 
where there are maximum Lindblad torques (~ H). This 
occurs when 

3 
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H 



(5) 



In this case, the usual torque cutoff condition for Lind- 
blad resonances is questionable. This argument suggests 
that the torque density maximum for Lindblad reso- 
nances should occur at a radial distance from the planet 



max (Rr, H) 



(6) 



When this condition is satisfied, the overall torque on 
the planet will be reduced, even if i?n < H, since reso- 
nances that lie closer than distance H from the planet are 
suppressed 3 . For the disk parameters considered in this 
section, this condition is satisfied for M p > 4 x 10 -4 M s 
(or 0.4 Mj). 

2.2.3. Saturation Effects of Coorbital Torques 

The flow in the coorbital region is trapped in horseshoe 
orbits. For a time-reversible system (e.g., no dissipation 
or migration), the streamlines are exactly periodic and no 
net torque occurs on the planet due to the disk (i.e., the 
torque saturates), except for possible initial transients 
due to initial conditions. However, turbulent viscosity 
introduces irreversibility that can lead to a net torque. 
The condition for saturation within the framework of the 

3 They may still partially contribute, due to their finite widths. 
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Fig. 1. — Torque per unit disk mass on the planet as a function 
of radius in units of the planet's semi-major axis, a. The vertical 
scale is in units of GM B (M P / M s ) 2 /a. The solid, long-dashed, dot- 
dashed, and short-dashed curves are for lMj;, 10 Me, 0.3 Mj, and 
1 Mj mass planets, respectively. The disk is modeled as a three- 
dimensional system. The vertical disk thickness is H/r = 0.05 for 
all the cases. Torque distributions are averaged over one orbital 
period. 

a-disk model is that the libration timescale of the fluid 
in the coorbital region is shorter than the viscous radial 
diffusion timescale across this region. Based on scaling 
arguments, the saturation condition is given by (Ward 
1992) 

<.*«""(£)"". P) 

For the parameters in this section, this constraint implies 
that for planets of order 10 Me or greater, the corotation 
torques should be saturated (small). Saturation effects 
should be important for the larger planet masses we con- 
sider. 

2.3. Numerical Results 
The torque per unit disk mass is defined by 
dT I 1 f 27r f°° \ 

where E(r, t) is the axisymmetric disk density (i.e., the 
surface density averaged over the azimuth (f>) and nota- 
tion (X(t)) is defined below equation (3). 

Numerically, the torque distribution per unit disk mass 
is determined by dividing the (three-dimensional) disk 
into a series of concentric shells, of radius R and thickness 
AR, centered at the origin and calculating the torque ex- 
erted by the shell and the mass of the shell. The torque 
per unit disk mass is obtain from the ratio of these two 
quantities 4 , averaged over an orbit period. We use the 
radial grid spacing on the base grid for the value of AR. 
The torques arising from within the Hill sphere of the 
planet are ignored in this section, but are included in 
later sections of this paper. We ignore such considera- 
tions here in order to compare results with the standard 
theory of coorbital and Lindblad torques, which does not 
include such contributions (Tanaka et al. 2002). 

4 There is a slight error of order (H/r) 2 in this procedure due 
to the difference between the spherical coordinate system used in 
the calculations and the cylindrical coordinates that apply to the 
definition of the torque in equation (3). 



The torque per unit disk mass for four planet mass 
cases is shown in Figure 1. The plots are normalized 
such that the torque densities in the four cases would 
be the same, according to linear theory, if the axisym- 
metric disk density gradients and gas properties (sound 
speeds and viscosities) were the same. That is, the torque 
density per unit disk mass is scaled by the square of 
the star-to-planet mass ratio. The 1 Me (solid line) and 
10A/e (long-dashed line) cases nearly exactly overlap as 
predicted, while the 0.3 Mj (dot-dashed line) and 1 Mj 
(short-dashed line) cases have a smaller scaled torque 
density. The scaling in the plot masks the fact that the 
results span a large range of parameter space. In going 
from 1 Me to 1 Mj there is a change in torque density 
by a large factor, 10 5 , while the discrepancy is about a 
factor of 2.5. 

The deviations in the 0.3 Mj and 1 Mj cases could be 
due to the modified torque cutoff, pressure gradients, and 
nonlinearities. Since i?H ^ H in these cases, Lindblad 
resonance contributions are weakened by the modified 
torque cutoff, as discussed in Section 2.2.2. Pressure 
gradients cause shifts in the resonance locations. For 
mild pressure gradients that change sign across the or- 
bit of the planet (as would occur for a mild gap), the 
resonances shift away from the orbit of the planet (see 
eq. 26 of Ward 1986). The shift would then cause the 
torques per unit disk mass to be weaker, as seen in the 
figure. The situation is more complicated in the case of 
stronger pressure gradients, as may occur for deep gaps, 
and the sign of the effect on the torque depends on the 
detailed shape of the density profile. Nonlinearities may 
play a role in the 1 Mj case, since there are shocks in the 
disk in that case, due to the strong forcing. But the to- 
tal torque is not expected to be substantially effected by 
nonlincarity. For a fixed smooth background disk den- 
sity distribution, resonant torques are quite insensitive 
to the level of nonlinearity (Yuan & Cassen 1994). For a 
1 Mj planet and a resonance with azimuthal wavenumber 
m = 20 = H / a, the nonlinearity is mild with nonlinear- 
ity parameter / = 0.6, as defined by Yuan & Cassen 
(1994). Some broadening of the torque density profile is 
predicted, while the total torque is reduced by only about 
1%. For much stronger nonlincarity, / = 3, the torque 
reduction is only 5%. This estimate is based on consid- 
ering only a single resonance. Many resonances overlap, 
increasing the level of nonlinearity. However, the theory 
does not describe overlapping resonances. So, although 
we cannot be definite about the importance of nonlinear- 
ities, indications for a single resonance suggest that they 
are not important. 

The torque density per unit disk mass for the 1 Mj 
planet in Figure 1 (short-dashed line) shows indications 
of saturation for \r — a\ < i?H- As discussed above, this 
effect is suggested by theoretical considerations. The 
torque density peak for the 1 Mj case is slightly displaced 
away from the planet relative to the smaller mass cases 
and lies close to a distance i?H — 0.07 a from the planet. 
This result is consistent with equation (6) in the 1 Mj 
case, \r - a\ ~ 0.07 a = 1.4 if. 

Figure 2 shows that the torque in the 1 Mj case is ac- 
quired close to the planet, well within the gap region. 
Most of the torque is accumulated by material with in- 
termediate/low density interacting with an intermediate 
magnitude torque per unit disk mass. About 80% of the 
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Fig. 2. — Azimuthally averaged surface density (long-dashed 
curve), disk torque per unit disk radius exerted on the planet 
(short-dashed curve), and cumulative torque (solid curve), i.e., 
torque per unit radius integrated outward, as a function of radius 
for a 1 Mj planet on a fixed circular orbit. The disk is modeled 
as a three-dimensional system. The unit of radius is the planet's 
orbital radius a. The surface density and cumulative torque are 
normalized by their absolute values at r = 2. The disk torque 
per unit disk mass is normalized by 10 3 GM s (M p /M s ) 2 /a. The 
plotted values are averaged over one orbital period. 

torque is due to material within a radial distance of 0.25 a 
from the planet. 

3. MIGRATING AND GROWING PLANETS 

We investigate the orbital migration of a planet that 
is undergoing run-away gas accretion. We consider sev- 
eral disk configurations, by changing the initial surface 
density, the pressure scale height, and the kinematic vis- 
cosity. We use disk models and numerical procedures 
similar to those introduced in section 2.1. Throughout 
this section, the disk is modeled as a three-dimensional 
system. The origin of the coordinate system is taken to 
be the star. The coordinate system rotates about the 
origin at a rate equal to the rotation rate of the planet 
around the star. We integrate the equations of motion 
of the planet, under the action of disk torques and ap- 
parent forces arising from the rotation of the reference 
frame, as described in DBL05. The unit of length is the 
initial star-planet separation ao (or 5.2 AU when convert- 
ing into physical units). The unit of time is the inverse 
of r2o, the initial angular speed of the planet. The unit 
of mass is the stellar mass M s (1 M©). 

The grid system achieves a linear base resolution of 
AR = ao AO = ao A(j> = 0.014 ao. In the coorbital re- 
gion around the planet, the linear resolution is about 
9 x 10~ 4 ao. Nested grid levels cover extended radial re- 
gions of the disk so that the planet remains within the 
domain covered by the most refined grid level over the en- 
tire orbital evolution. Convergence tests were carried out 
with a grid system that used a volume resolution (3/2) 3 
times as high throughout the whole disk domain and on 
all grid levels. No significant differences are observed 
(see Appendix A.l). To avoid depletion of the disk inte- 
rior of the planet's orbit, we apply nonreflecting bound- 
ary conditions to the inner grid (radial) border. We test 
our results against possible boundary condition effects in 
Appendix A. 2 by applying outflow boundary conditions 
and moving radial disk boundaries farther away from the 



planet's orbit in both directions. No important effects are 
observed. Near the planet we apply accreting boundary 
conditions on the g described in section 3.1.1. We 
consider planetary mass increases that extend over more 
than two orders of magnitude and a range of disk surface 
densities. 

To avoid possible spurious torques exerted by material 
gravitationally bound to the planet, contributions from 
within Rh/2 of the planet are not taken into account. We 
report in Appendix A. 3 on the sensitivity of the results to 
the radius of the excluded region by considering a smaller 
radius. We find that the changes are not significant. 

We generally initiate the calculations with a planet 
mass M p = 1.5 x 10~ 5 M S , or 5M E . However, in some 
applications discussed in section 4, we use an initial mass 
M p = 3 x 10~ 4 M s (about 0.3 Mj) in order to study the 
effects on migration of releasing a more massive planet 
in an unperturbed disk. 

3.1. Planet Mass Growth 
3.1.1. Gas Accretion 

In the core accretion scenario of giant planet forma- 
tion, prior to the phase of run-away gas accretion, the 
rate at which gas is accreted is largely determined by 
the ability of a planetary core's envelope to radiate away 
the energy delivered by gas and solids (phase of slow 
gas accretion, see e.g., Hubickyj et al. 2005). During 
the initial stages of planet growth, the accretion of solids 
dominates, and the dissipation of the kinetic energy of 
the impacting solids provides an important heat source 
for the accreted gaseous envelope. Models of Hubickyj 
et al. (2005), which ignore the effects of planet migra- 
tion, experience a depletion of solid disk material in the 
vicinity of the planet and consequently a reduction in 
the envelope heating rate. When the mass of the gas (in 
the envelope) is comparable to the mass of solids (in the 
core), the pressure gradient cannot prevent the gravita- 
tional collapse of the envelope. This situation results in 
a sudden increase of the gas accretion rate and a rapid 
growth of the planet's mass, the so-called run-away gas 
accretion phase (e.g., Wuchterl 1993; Pollack et al. 1996). 

The models presented here assume run-away gas ac- 
cretion. They do not account for the thermal struc- 
ture and detailed microphysics of a planet's envelope. 
Therefore, we do not determine self-consistent gas accre- 
tion rates, prior to the phase of run-away gas accretion 
(Mp < 10 .Me). The models also ignore the effects of 
heating by impacting solids that act to slow the gas ac- 
cretion, as the planet migrates out of the region of de- 
pleted solids. During the run-away gas accretion phase, 
the accretion rate onto the planet is only limited by the 
amount of gas that the disk is able to supply. The calcu- 
lations described here provide estimates of such limiting 
gas accretion rates during the run-away gas accretion 
phase. 

In these models, we adopt a prescription that gas 
within a distance of i? a cc = 0.1 i?n from the planet can 
accrete onto it. Accreted gas is removed from the disk 
and its mass is added to the planet mass. For the mod- 
els we consider, this distance is safely smaller than the 
possible characteristic accretion radii: the Hill radius, 
-Rh, arid the Bondi radius, Rb (distance beyond which 
the thermal energy of the gas is larger than the gravi- 
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Fig. 3. — Mass evolution of a protoplanet having initial 
planet mass 5 undergoing run-away gas accretion in a three- 
dimensional disk with initial surface density E p = 3 X 10 -4 M a a^ 2 
or about 100 gem -2 at the planet's initial orbital radius of 5.2 AU 
(solid line) and S p = 9 X 10 — 4 M S a^" 2 or about 300 gem -2 (dashed 
line). In both cases, the disk thickness is H/r = 0.05 and the tur- 



bulent viscosity parameter is u = 1 X 10" 
5.2 AU). The time refers to orbits at ao 
years. 



a 2 Q (a = 0.004 at 
5.2 AU or about 12 



tational energy that binds the gas to the planet). The 
distance i? acc is at least a factor of 3 smaller than Rb- 
Therefore, this mass removal prescription should not de- 
termine the accretion rate for the case of run-away gas 
accretion (see also Tanigawa & Watanabc 2002). The 
amount of material accreted per time-step At is given 
by (Ai/r acc ) / pdV, where dV is the volume element and 
r acc is a removal timescale. The integral is performed 
over the sphere of radius 0.1 Rb centered on the planet. 
Here we set r acc = 0.1 fl^ 1 within the sphere of radius 
0.05 i? H and r acc = 0.3 il^ 1 for 0.05 Rb < S < 0.1 Rn (S 
is the distance from the planet). 

3.1.2. Mass Evolution 

In this section we describe the accretion rates of mi- 
grating, mass-gaining planets. Figure 3 shows the planet 
mass as a function of time, M p = M p (t), for a model 
with initial (unperturbed) surface density at the initial 
orbital radius of the planet E p = 3 x 10~ 4 M s a^ 2 [solid 
line). For a planet orbiting a Solar mass star at 5.2 AU, 
this density is about 100 gem -2 , roughly corresponding 
to the minimum mass solar nebula. 

The mass evolution can be understood in terms of 
Bondi and Hill accretion. Consider a simple model in 
which gas is captured within some radius, S c , of a planet 
and assume S c < H . Mass is accreted with some velocity 
relative to the planet of order Q S c , and so the mass ac- 
cretion rate in a three-dimensional disk (where p ~ E / H) 
is estimated as 

M p ~^QS 3 c , (9) 

where we take S c as either the Bondi or Hill radius, with 
the Bondi radius given by Rb = GM p /c 2 s and the Hill 
radius given by i? H = a [Af p /(3 A^)] 1 / 3 . 

In the case that gas pressure prevents the gas from be- 
ing bound to the planet within the Hill sphere (or, equiv- 
alently, that pressure forces dominate over gravitational 
three-body forces), we expect the Bondi description to 



Fig. 4. — Mass growth rate 1/tq = M p /M p in units of inverse 
orbital periods at the initial radius of the planet, Qo/(2n), plotted 
against M P /M B for the solid curve case in Figure 3. The dashed 
line plots the growth rate according to equation (15). The slopes 
of the two dashed line segments are predicted by the model. The 
two free parameters, Cb and Ch, are dimensionless constants of 
order unity that control the intercepts and are fit to the solid curve. 
The slanted portion of dashed line corresponds to accretion within 
the Bondi radius, given by 1/tb in equation (15) with Cb = 2.6. 
The horizontal portion of the dashed line corresponds to accretion 
within the Hill radius for a disk with no gap, given by 1/tjj in 
equation (15) with Ch = 0.89. At higher planet masses, the growth 
rates drop due to the presence of the tidally produced gap. 



be appropriate. This condition is that 

> GM P 



Rh 



or 



(10) 
(11) 



Therefore in the general case we take 

5 c = min(i? B ,i? H ). (12) 

It then follows that the Bondi and Hill mass growth 
rates, M p /M p , of the planet are given by 



l/r B = C B ft 



So 2 

i s« 2 



(I) 7 
(I)- 



Mp 



(13) 



(14) 



where Cb and Ch are dimensionless coefficients of order 
unity. The overall mass growth rate is given by 



where 



1/tg 



M, 



1/t b for M p < A/ t 
1/th for M p > M t 




(15) 



(16) 



is the transition planet mass where th = Ta- 
in Figure 4, we plot the mass growth rate, 1/tg, for 
the solid curve case in Figure 3. We applied equa- 
tion (15) and adopted constant values of E = E(ao), 
at time t = 0, and r2 = Slo- The figure shows that 
the Bondi and Hill accretion rates in equation (15) agree 
with the simulation results for values of Cb = 2.6 and 
Ch = 0.89. The transition mass in this case evaluates 
to Af t = 4.2 x 10~ 5 A/ S . It lies between the Bondi and 
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Hill accretion regimes in the figure, at the intersection 
between the two dashed line segments. For larger values 
of planet mass, M p > 2 x 1CT 4 M S « 4.8 M t , this simple 
estimate of the mass growth rate breaks down because 
the density is depleted near the planet due to the onset 
of gap formation. The density near the planet is reduced 
by about 40% when M p w 2 x 1CT 4 M s (see Fig. 6, right 
panel). In addition, the Hill radius becomes comparable 
to H , since i? H = H = 0.05 a for M p = 3.75 x 10" 4 M s . 

Simulations carried out in two dimensions would have 
different scaling behavior, since the right-hand side of 
equation (9) would be EOS' 2 . The dependence of the 
mass accretion rate on planet mass and disk sound speed 
then artificially deviates from the three-dimensional case. 
In two dimensions we have that 1/tb oc (M p /M s ) (a/H) 
and 1/th oc (M s /M p )^ 3 . 

The maximum of the accretion rate for the solid curve 
case of Figure 3 is M v ~ 5 X 10~ 3 £ p a 2 ~ 1.5 x 10" 3 Mj 
per orbit and occurs when M p w 0.3 Mj. This re- 
sult is consistent with the previous findings of D'Angelo 
et al. (2003) and Bate et al. (2003), who considered 
planets on fixed orbits. Also displayed in Figure 3 is 
the planet's mass evolution in a disk with initial S p = 
9 x 10 _4 M s aQ 2 [dashed line) or about 300gcm~ 2 at 
5.2 AU. For M p /M s < 10 -4 , the accretion rate is a fac- 
tor of 3 larger than that of the lower density disk case 
(solid line). Hence, equation (15) applies to the growth 
rate with the same coefficients Cb and Ch as those given 
above. For larger planet masses, the accretion rate keeps 
increasing until M p sa 0.7 Mj, at which point M p starts 
to decline very rapidly as M p grows further. This is 
because effects due to gap formation are delayed. The 
timescale required to form a gap of half-width £i?H is 
7"ga P ~ £ 5 q -1 / 3 (see, e.g., Bryden et al. 1999), where 
£ ~ 2 (see long-dashed line in Fig. 2). In the lower 
density disk model (solid curve in Fig. 3), r gap < tq 
for M p /M s > 10~ 4 . In the higher density disk model 
(dashed curve), T gap becomes shorter than tq only when 
M p > 0.7 Mj. 

In Figure 5, the mass evolution is shown for cases 
in which E p = 3 x 10" 4 M s a 2 w 100 gem" 2 , but 
with different scale heights, H, and kinematic viscosi- 
ties, v. Near M p = lMj, the accretion rates of the 
two models with different H/r (solid and long-dashed 
lines), but the same S p and v, are nearly equal, with 
M p w 3 x 10~ 3 S p a 2 ~ 9 x 10~ 4 Mj per orbit. At larger 
planet masses, M p is smaller in the case of a colder disk 
(long- dashed line) because of the stronger tidal torques 
exerted by the planet on the disk material that produce 
a wider gap. When M p w lMj, the simulation with 10 
times larger viscosity (short- dashed line) yields an ac- 
cretion rate that is a factor of nearly 8 larger. This re- 
sult is consistent with previous two-dimensional studies 
of planets on fixed orbits that do not gain mass. For 
M p w 1 Mj these studies showed that M p scales ap- 
proximately linearly with vTi, the overall disk accretion 
rate evaluated just outside the gap (Klcy 1999; Lubow & 
D'Angelo 2006). 

3.1.3. Mass Within the Hill Sphere 

We discuss here the relevance of torques exerted on a 
planet and originating within the planet's Hill sphere. 
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Fig. 5. — Mass evolution of a protoplanet having initial 
planet mass 5 Me and undergoing run-away gas accretion in a 
three-dimensional disk with initial surface density E p = 3 X 
U)- 4 M 8 a, 







100 gem -2 at the planet's initial orbital radius 
ao = 5.2 AU. The solid line represents a case with H/r = 0.05 
and turbulent viscosity i/=ix 10~ 5 a 2 Qo (ot = 0.004 at 5.2 AU), 
the long-dashed line refers to a variant model with H/r = 0.04 
and the same v value (a = 0.006 at 5.2 AU), and the short-dashed 
line represents a variant model with v = 1 X 10~ 4 a 2 , f!o (a = 0.04 
at 5.2 AU). The time refers to orbits at ao = 5.2 AU or about 12 
years. 

We may expect that material gravitationally bound to 
the planet should not be capable of exerting significantly 
strong torques, if resolution is appropriate (DBL05). In 
some situations, if the local density is large, any torque 
imbalance can be easily amplified by lack of numerical 
resolution (because torques depend on 1/S 2 , where S is 
the distance to the planet). Artificial effects may arise 
when the mass within ~ i?n of the planet is larger than 
the planet's mass. However, not all this material is nec- 
essarily bound to the planet. Because of the nonspherical 
nature of the Roche lobe, the Hill radius represents an 
overestimate for the size of the region where gas is bound 
to the planet (Paczyhski 1971; Egglcton 1983). We have 
found that accumulated gas may be bound to the planet 
within distances shorter than -Rh/2 from the planet (see 
Appendix D.l). 

In all the cases discussed in this section, the amount 
of material that lies within -Rh/2 of the planet is smaller 
than M p , throughout the evolution, by several orders 
of magnitude. For models in Figure 3, as well as for 
those in Figure 5, the ratio of these two masses ranges 
from less than ~ 10 -3 to ~ 10~ 2 , depending mainly on 
the planet's mass. We also consider models with initial 
densities larger than those discussed here (described in 
section 4). However, this mass ratio remains on the or- 
der of 10~ 2 or smaller. Therefore, due to the accretion 
boundary condition employed here at the planet location, 
these models do not experience a build-up of mass near 
the planet (with possible effects on planet migration). 
The accreted mass is accounted for by the increase in 
the planet mass. 

3.2. Planet Migration 

3.2.1. Theoretical Regimes of Migration 

A planet that grows in mass from a few Earth-masses 
to a few Jupiter-masses is susceptible to two "classical" 
regimes of migration. The Type I regime is expected 
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when the planet causes small, linear disk density pertur- 
bations (e.g., Ward 1997; Tanaka et al. 2002). In the op- 
posite limit, Type II occurs when the planet mass is large 
enough to cause nonlinear density perturbations that re- 
sult in a density gap along its orbit (Lin & Papaloizou 
1986). 

For the parameters we adopt (pressure scale height 
H/r ~ 0.05, kinematic viscosity of disk v > 1 x 
10~ 5 a§ fi 0j an d initial planet mass M p /M s = 1.5 x 10~ 5 
(or M p = 5 Me), it is expected that the initial evolu- 
tion of the planet will follow Type I migration, since the 
usual gap opening criteria are not satisfied. In the lin- 
ear theory of Tanaka et al. (2002), the rate of migration 
resulting from the action of both Lindblad and (unsatu- 
rated) coorbital corotation torques is given by 

^ = -( 2 - 73 + L08s K^) ^ (17) 

where s is the slope of the unperturbed surface den- 
sity. For the case of saturated (zero) coorbital corotation 
torques, the migration rate is given by 

* = -<« 8 -o,o»> (tf)> s <v m 

The conditions for saturation are discussed in sec- 
tion 2.2.3. For higher planet masses that arise in the 
later stages of the simulations, the torques are expected 
to be saturated. 

In the presence of a sufficiently clean density gap and 
for a planet whose mass is less than the local disk mass, 
the rate of migration follows Type II theory that is dic- 
tated by disk viscous inflow 



Note that if there is residual material in the horseshoe 
orbit region, the migration rate can differ from that in 
equation (19). The coefficient £ on the right-hand side of 
equation (19) is of order unity and also depends on the 
evolutionary state of the disk. For a steady-state disk, 
the coefficient is 3/2. But for nonsteady disks where vY, 
varies in radius, as in our initial states, the coefficient 
may differ by order unity amounts. 

In the unsaturated case, some nonlinear effects of the 
corotation resonance can cause migration rates to dif- 
fer from those predicted by equation (17) (Masset et al. 
2006). For s = 1/2, H/r = 0.05, these effects occur in 
the range of masses is between w 10 Me and « 20A/e- 
However, in the models presented here, the planet grows 
too quickly through this mass range (taking less than a 
few tens of orbits) to significantly affect migration (see 
Fig. 23 in Appendix B). 

When the amount of material in the horseshoe orbit 
region is larger than the planet's mass, a regime of fast 
migration known as Type III may occur. The origins 
of such a regime are not yet entirely clear. The model 
of MP03 suggests that it is driven by strong corotation 
torques originating from material that streams past the 
planet, while the planet is moving in the radial direction. 
However, an analytic model of OL06 suggests that such 
torques could originate from trapped librating gas. A 
somewhat similar model was developed by Artymowicz 
(2004). 



3.2.2. Orbital Radius Evolution 

We evaluate quantities S p , H, and f2 p at the planet's 
orbital radius, a. Surface density S p = T, p (a) is eval- 
uated according to its initial value T, p (a) oc (ao/a) s , 
and so ignores evolutionary effects and tidal gap forma- 
tion. The planet mass M p is regarded as a function of 
time that we obtain from our simulations, via piecewise 
polynomial fits. For the numerical models we consider, 
s = — dlnEp/dlna = 1/2. Equations (17) and (18) are 
then solved numerically, providing the migration tracks 
ai = ai(t). 

In the left panel of Figure 6, we compare such tracks 
with outcomes from our simulations. For the first 400 
orbits, while i? H < 0.9 H and M p < 0.27 Mj, the orbital 
radius (i.e., semi- major axis) evolution is in good agree- 
ment with the results of Type I migration. The unsatu- 
rated coorbital torques appear to give a better fit than 
the saturated ones. But this is not always the case, as we 
see later when different disk parameters are considered. 
The right panel of Figure 6 plots the density evolution 
of the gas near the planet, S^, computed as ratio of the 
disk mass in the radial band \r — a\/a < H/r to the area 
of the band (Y.% is the local initial value of £b). It shows 
that the migration rate follows the Type I tracks on the 
left while the disk density near the planet remains close 
to the local initial disk value, assumed in equations (17) 
and (18). Up to a time of about 400 orbits, the density 
near the planet is reduced below its local initial value by 
less than 20%. At time of about 600 orbits, the density 
near the planet's orbit is reduced by about a factor of 3, 
and we should expect the migration rates deduced from 
the simulation to be substantially slowed below the rates 
based on Type I theory, in accord with the results on the 
left panel. After about 1000 orbits, when M p > 0.9 Mj, 
the migration rate in the simulation becomes comparable 
to the (local) viscous inflow rate (long-dashed line). At 
this point, the disk density near the planet is depleted 
by a factor of about 30. 

The torque per unit disk mass as a function of dis- 
tance from the planet for the case in Figure 6 is plotted 
in Figure 7. The plot shows very similar behavior to the 
case of a stationary, nongrowing planet seen in Figure 1. 
Therefore, there is no evidence that planet migration or 
growth substantially affects the disk-planet torques for 
these model parameters. In particular, there is no evi- 
dence for strong coorbital torques. 

The results obtained from a model with H/r = 0.04 
(i.e., with a lower disk temperature compared to the 
model in Fig. 6) are shown in the left panel of Figure 8. 
As in the case of the warmer disk, the Type I migra- 
tion tracks (short-dashed curves) reproduce reasonably 
well the radial migration from the simulation (solid line) 
while M p < 0.14Mj (see long-dashed line in Fig. 5) or 
Rk ;$ 0.9 H. As before, the right panel of Figure 8 shows 
that the migration rate follows the Type I tracks while 
the disk density near the planet remains close to the 
local unperturbed value. Again, when M p > 0.75 Mj, 
Eb/E'q < 0.03 and \da/dt\ is on the order of the viscous 
inflow velocity (long-dashed line). 

The dependence of migration on viscosity was investi- 
gated by running a simulation with kinematic viscosity 
v = 1 X 10~ 4 aQf2o (ol = 0.04), ten times the value in 
Figure 6 with all other parameters being the same. The 
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Fig. 6. — Orbital migration of a planet undergoing run-away gas accretion. Left: Orbital radius in units of ao (5.2 AU), as a function of 
time in units of the initial orbital period (fa 12 years). The initial planet mass is 5 Mj;. The initial surface density is S p = 3 X 10 -4 M a aZ 2 Rj 
100 gem -2 at the planet's initial orbital radius and H/r = 0.05. Solid curve: Results from the three-dimensional numerical simulation of 
a migrating, gas-accreting planet. Short-dashed curves: Predictions based on Type I migration theory, obtained by solving equations (17) 
and (18), for a planet that undergoes the mass growth given by the solid line in Figure 3 and is embedded in a disk with the initial 
unperturbed density distribution. The upper (lower) curve is for migration with unsaturated (saturated) coorbital torques. Long-dashed 
line: Consistent with Type II migration, the line has slope — 1.5f/a and passes through a ss 0.8ao when M p as 0.9Mj. Right: Average 
disk density near the planet relative to the local initial value as a function of time. The density is averaged over a band of radial width 
2H centered on the orbit of the planet (see text for details). Solid circles mark times when the mass ratio M p /M s is equal to 5 X 10 — 5 
(M p = 16.7 M B ) and when it is an integer multiple of 1 X 10~ 4 {M p = 33.3 M E ). 
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Fig. 7. — Torque per unit disk mass on the planet as a function of 
normalized distance from the migrating and growing planet plotted 
in Figure 6 (S p = 3 X 10~ 4 M s ag 2 M 100 gem -2 at the planet's 
initial orbital radius and H/r = 0.05). The vertical scale is in units 
of GM s (Mp/M s ) 2 /a, where a = a(t). The solid, long-dashed, dot- 
dashed, and short-dashed curves refer to times when M p = 6.0 Me, 
9.3 M E , 0.36Mj, and 1.0 Mj, respectively. 

results are shown in Figure 9. The left panel shows the 
orbital migration from the simulation as a solid curve and 
the Type I migration based on equations (17) and (18) 
as dashed curves. In this case, the relation M p = M p (t) 
represented by a short-dashed line in Figure 5 is used 
in equations (17) and (18). The long-dashed line indi- 
cates a migration at a constant rate of \a\ ~ 0.7 v/a, 
with a ~ 0.9 ao. The long-dashed line passes through a 
range of masses that spans from w 0.2 Mj to « 1.2 Mj. 
However, at M p sa lMj (t as 600 orbits), the density 
gap along the planet's orbit has not yet fully formed. 
This can be observed on the right panel of Figure 9, 
which displays the averaged disk density near the planet 
normalized to the local unperturbed (initial) disk value. 



There is a drop of only a factor of 2.5 in the disk den- 
sity near the planet by the time M p as 1 Mj. The reason 
is that one of the conditions for steady-state gap forma- 
tion, M p /M s > 40is/(a 2 n) ~ 4x 10~ 3 (Lin & Papaloizou 
1993), is not fulfilled in this higher viscosity case until 
M p > 4Mj. At about 780 orbits, Eb/E'b ~ 0.1 but 
the planet mass has reached beyond 2 Mj and is there- 
fore more massive than the local disk mass. At those 
stages of the orbital evolution, inertia effects and further 
gap clearing are likely playing an important role in re- 
ducing the migration rate, as demonstrated in the next 
paragraph. 

Figure 10 displays a comparison of the orbital radius 
evolution from two calculations. The solid line is the 
same as that in the left panel of Figure 9. The dot- 
ted line with solid circles is the outcome of a three- 
dimensional simulation in which the planet mass is fixed 
at M p = 1 Mj . Material is removed from the vicinity 
of the planet according to usual the procedure we apply 
(see section 3.1.1), but in this case it is not added to the 
mass of the planet. The planet's orbit is held fixed for the 
first 100 orbital periods, after which time it is allowed to 
evolve under the action of disk torques. The plot shows 
that there is general agreement, while M p ~ 1 Mj, with 
the variable mass model and that the effect of adding 
mass to the planet in this regime is to slow its migration 
rate. 

The local viscous timcscalc, t v = r 2 /v, in the models 
presented in Figures 9 and 10 is about 1600 orbital peri- 
ods at r = ao. Therefore, one might wonder whether the 
viscous evolution of the disk at radii larger than the outer 
grid boundary has any significant impact on the orbital 
evolution of the planet. We address this issue in Ap- 
pendix B and show that extending the disk further out at 
larger radii does not affect the migration tracks shown in 
Figures 9 and 10. In Appendix B, we also present results 
for cases with viscosity parameter a = 0.2 (kinematic 
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Fig. 8. — Orbital migration of a planet undergoing run-away gas accretion. Same as Figure 6, but for a cooler disk with aspect ratio 
H/r = 0.04 (the same v and initial E p ). Left: The theoretical Type I migration tracks (dashed curves) use the mass evolution shown as a 
long-dashed curve in Figure 5. As in Figure 6, the upper (lower) short-dashed curve is for unsaturated (saturated) coorbital torques. The 
long-dashed line, representing Type II migration, has a slope equal to — 1.5 v/a and passes through a » 0.85 ao, when M p m 0.9 Mj. Bight: 
Normalized disk density near the planet as a function of time, as described on right panel of Figure 6 (see also text). Solid circles mark 
times when M p /M s is 5 X 10~ 5 (M p = 16.7 M E ) or an integer multiple of 1 X 10~ 4 (M P = 33.3 M E ). 



viscosity v = 5 x 10~ 4 ag^o) that have t v ~ 320 orbits 
at r = ao . This case also leads to inward migration that 
can be interpreted as a Type I regime, partially modified 
by the perturbed surface density of the disk. 

4. TYPE III MIGRATION 

Figures 6 and 8 indicate that a growing planet un- 
dergoes Type I migration, as long the disk density near 
the planet remains undepleted. At higher planet masses 
where the gap opening sets in, there is a smooth tran- 
sition towards Type II migration with migration speeds 
that are on the order of the viscous inflow velocity. There 
is no evidence for another form of migration, since the 
torque distributions are essentially the same in the mi- 
grating and nonmigrating cases explored thus far (com- 
pare Figures 1 and 7). Type III migration was suggested 
to involve coorbital material that provides a fast form 
of migration (MP03). In this section we discuss planet 
migration for several variants on the models of section 3 
that should be favorable for a Type III regime of migra- 
tion. We describe a case that appears to exhibit Type III 
migration. 

4.1. Higher Disk Mass 

Coorbital torques are stronger for higher mass disks. 
Masses in the coorbital region arc on the order of 
8ir i?H a S(a). For the model presented in Figure 6, 
involving disks of relatively low density, the coorbital 
disk mass is approximately equal to the planet mass 
when M p w 0.2 Mj. We describe here results of three- 
dimensional calculations with initial surface densities 
S p = 9 x 10~ 4 A/ S ao 2 w 300gcm" 2 and S p = 1.5 x 
10~ 3 M s a^ 2 w 500 gem -2 at the planet's initial orbital 
radius of ao = 5.2 AU. The mass evolution in the for- 
mer case is plotted as the dashed line in Figure 3. The 
mass evolution in latter case is similar, but the growth 
proceeds very rapidly reaching about 1 Mj within 130 
orbital periods. The resulting orbital radius evolution 
for both simulations is plotted in Figure 11 (left panel) 
along with the average disk density near the planet nor- 
malized to the local unperturbed value (right panel). 



For both cases presented in the figure, at earlier times 
(t < 170 and t < 100 initial orbits, respectively), the 
simulated migration rates are comparable to the Type I 
rates. During that stage of the evolution, the coorbital 
region is more massive than the planet. In the model 
with initial S p « 300gcm~ 2 at 5.2 AU (upper migration 
track in Fig. 11), for times t < 170 orbits (M p < 0.3 Mj) 
the coorbital region mass to planet mass ratio is larger 
than 2. In the model with initial S p as 500gcm~ 2 
at 5.2 AU (lower migration track Fig. 11), for times 
t < 100 orbits (M p < 0.4 Mj) the ratio of coorbital region 
mass to planet mass is larger than 3. However, during 
those stages, the results are generally consistent with the 
Type I migration and some slowing at later times, with 
no indication of another form of migration. 

To examine the situation in more detail, we plot the 
torque per unit disk mass as a function of distance from 
the planet in Figure 12. The plot shows very similar be- 
havior to the case of a nonmigrating, nongrowing planet 
seen in Figure 1, as well as to the case of a migrating, 
growing planet within a lower density disk presented in 
Figure 7. Again, there is no evidence that planet mi- 
gration or growth substantially affects the disk-planet 
torques for the parameters adopted in these models. 
Furthermore, there is no evidence for strong coorbital 
torques dominating planet's migration. 

In carrying out calculations at higher disk masses, we 
have introduced a possible inconsistency between the or- 
bital motion of the disk and the planet. The orbital mo- 
tion of the planet is affected by the axisymmctric grav- 
itational force of the disk. On the other hand, the mo- 
tion of the disk near the planet is not affected by this 
force, since disk self-gravity is ignored. This difference 
in rotation rates can lead to an artificial increase in the 
planet migration rate (Pierens & Hure 2005; Baruteau & 
Massct 2008). This issue has some quantitative effect on 
our results in this section. But, the qualitative results 
(approximately following the expectations of standard 
Type I and II theory) remain. We examine this issue 
further in Appendix C. 
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Fig. 9. — Left: Same as left panel of Figure 6, but for a disk with ten times the turbulent kinematic viscosity (v = 1 X 10 -4 a 2 , Qo or 
a = 0.04), same H/r, and initial E p . As in Figure 6, the upper (lower) short-dashed line is for unsaturated (saturated) coorbital torques, 
using the mass evolution shown as a short-dashed curve in Figure 5. The long-dashed line representing Type II migration has a slope 
equal to —0.7v/a. Right: Average disk density near the planet relative to the local initial (unperturbed) value as a function of time, as 
in the right panel of Figure 6. Solid circles mark times when M p /M a is 5 X 10 — 5 (M p = 16.7Me) or an integer multiple of 2 X 10~ 4 
(Mp = 66.6 M B ). 
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Fig. 10. — Comparison of radial migration obtained from the sim- 
ulation on the left panel of Figure 9 (solid line) with that obtained 
from a similar three-dimensional simulation (dotted line with solid 
circles) with a fixed mass planet M p = 1 Mj (see text for further 
details). 

4.2. Higher Initial Planet Mass 

We have shown that if a low mass protoplanet is al- 
lowed to rapidly grow in mass while it migrates, the or- 
bital radius evolution begins at the Type I rate (eqs. 17 
and 18) and approaches the Type II migration rate as 
a clean gap develops. Since the evolving planet gains 
mass at the fastest possible rate, the run-away accretion 
rate, the time available for gap clearing is relatively short. 
Such conditions should be favorable for migration domi- 
nated by coorbital torques. But as we saw in Figure 11, 
such situations only reveal Types I and II migration. In 
this section, we explore a more extreme situation for pro- 
viding coorbital material. We consider the case that a 
planet of higher initial mass (higher than the 5 Me con- 
sidered thus far) is suddenly immersed in a smooth disk. 
Gap clearing is then not initially present for the higher 
mass planets. More coorbital gas is available for affecting 
migration. 

We consider a planet with initial mass M p = 0.3 Mj 



(M p /M s = 3 x 10 4 ) that is allowed to grow and migrate 
in a three-dimensional disk with initial density T, p 
300 g cm -2 at ao = 5.2 AU (same as the lower initial den- 
sity disk in Fig. 11, H/r = 0.05, and u=lx 10" 5 a§ fi )- 
Its orbital radius evolution is plotted as a dotted curve 
with solid circles in Figure 13, together with the migra- 
tion track of the model that starts with M p = 5 Me at 
t = (plotted as a solid curve). For purposes of com- 
parison, the initial orbital radius ao for the dotted curve 
case is chosen to be the a value of the solid curve case 
when its planet mass is also 0.3 Mj. Given the large 
initial mass of the planet for the dotted curve case, the 
mass growth is very rapid: the planet gains about 0.7 Mj 
over the first ~ 50 orbits of evolution. The figure shows 
that, with these disk conditions, migration rates differ 
only for a brief period of time, but they soon converge 
to values compatible with orbital migration in the more 
relaxed disk (compare slopes of solid and dotted lines). 
The dashed lines in the figure show that the planet ini- 
tially migrates at the Type I rate 5 . But it later slows 
to nearly the same rate as the solid curve case. Again, 
there is no indication of Type III migration. 

4.3. Nongrowing Planets in a Colder Disk 

We consider the case of a nongrowing 0.3 Mj planet by 
removing gas mass near the planet without adding the 
mass of this material to the planet's mass. This situation 
may mimic the effects of an efficient disk wind. These 
models differ from those in MP03 and DBL05, who con- 
sidered nonaccreting planets, only with respect to the 
accretion boundary conditions near the planet and the 
time of planet release. 

Unlike the mass removal case, the nonaccreting case 
may introduce a complication because of the buildup 
of gas within the planet's Hill sphere, which can be- 
come more massive than the planet. It has been ar- 
gued that inertia effects from material close to the planet 
could introduce complications in self-consistently analyz- 

5 Note that, for a constant mass planet and S p cx a -1 / 2 , it 
follows that | ai| on a (see eq. 17). 
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Fig. 11. — Right: Orbital evolution under the same conditions as the model in Figure 6, but with higher disk densities. Solid curves: 
Simulation results for orbital migration of a planet in a three-dimensional disk with initial surface density equal to S p = 9 X 10~ 4 M s a^ 2 , 
or about 300gcm~ 2 at ao = 5.2 AU (upper migration track), and S p = 1.5 X 10~ 3 M s a^ 2 , or about 500gcm~ 2 (lower migration track). 
Dashed curves: Predicted orbital migration according to Type I theory, equations (17) (upper curve of pair for unsaturated coorbital 
torques) and (18) (lower curve of pair for saturated coorbital torques). Flight: Average disk density near the planet relative to the local 
initial (unperturbed) value as a function of time, as in the right panel of Figure 6. Solid circles mark times when M p /M s is 5 X 10~ 5 
(M p = 16.7 Mg) or an integer multiple of 2 X 1CT 4 (M p = 66.6 Mg). 



ing the dynamics of the system (Papaloizou ct al. 2007) . 
Appendix D describes some effects of the nonaccreting 
boundary condition. To avoid this potential problem, 
we remove gas near the planet and ignore torques ex- 
erted by the gas on the planet within the inner half of 
the Hill sphere (by radius), where most of the bound gas 
resides (see Fig. 29). 

We are interested in seeing whether these situations 
could give rise to strong torques outside the Hill sphere 
that cannot be accounted for by Type I or Type II theory 
in the case of a planet of fixed mass. As we show below, 
there are conditions under which such strong torques oc- 
cur in the coorbital region. 

4.3.1. Simulations Setup 

We consider a Saturn mass (0.3 Mj) planet in a disk 
having H/r = 0.03. The initial (unperturbed) disk 
surface density varies as X = S p (a )(a /r) 3 / 2 , with 
^p(ao) = 2x 10 _3 A/ S a^ 2 « 670gcm -2 , and kinematic 
viscosity v = 1 X 10~ 5 af t Oo- Most of these calculations 
are carried out in two dimensions since i?H — 1.5 i?. 
However, we checked that results from three-dimensional 
models are in general agreement with those from two- 
dimensional models (see dotted curves in Fig. 14). Simu- 
lations in three dimensions use the grid system outlined 
in section 3. As in the three-dimensional case, the two- 
dimensional grid has a linear base resolution of 0.014 ao. 
In the coorbital region around the planet, the linear res- 
olution is 0.02 i?H- Since we intend to study some global 
properties of flow dynamics in the coorbital region, three 
grid levels extend 2tt in azimuth around the star. Con- 
vergence tests at these grid resolutions are presented in 
Appendix D. 

As anticipated above, here we assume that some pro- 
cess removes gas from the disk, according to the proce- 
dure detailed in section 3.1.1, but that the planet mass 
remains constant. The migration rate of the fixed mass 
planet is not substantially affected by the assumption 
that the gas is removed. In Appendix D we show that 



configurations with a nonaccreting planet result in simi- 
lar migration tracks. Hence, our conclusions would apply 
to nonaccreting planets as well. 

In MP03 and DBL05, the planet's orbital radius was 
initially fixed for over 470 orbits, so that a time-steady 
disk gap would form before it underwent migration. 
Here we reconsider that configuration, but examine cases 
where the planet's initial orbital radius is fixed for a 
shorter time, only 100 orbits. This case is somewhat like 
that of section 4.2 which has more gas in the coorbital 
region (lower curve case in Fig. 11), but instead has a 
fixed mass planet in a cooler disk. The following factors 
applied here should help increase the torques from the 
coorbital region: cooler disk, fixed planet mass, higher 
disk density, and reduced time on initially fixed orbit. 

4.3.2. Results 

Figure 14 shows that the migration timescale of the 
planet is quite short and that it lengthens as the re- 
lease time increases (and the gap deepens). We focus 
on the case with t r i s = 100 (solid line case in the figure), 
which has a migration timescale of order 100 initial or- 
bital periods. Although short, this migration timescale 
is longer than the Type I migration timescale that would 
be predicted if the planet did not open a gap (lower 
short- dashed curve of pair 6 ). The planet does open a 
partial gap, as seen in Figure 15. So it might appear 
that a weakened form of Type I migration, due to par- 
tial gap opening, could explain the simulated migration 
rate. However, we demonstrate below that the migration 
cannot be explained by the usual Type I theory. 

Figure 16 shows the torque distribution per unit disk 
mass exerted on the planet, as defined by equation (8). 
The torque density distribution at the time of release 
of the planet, 100 orbits after the start of the simula- 
tion, reveals a curve characteristic of Type I torques. 

6 Note that, for a constant mass planet and S p <x a -3 / 2 , it 
follows that a\ is a constant (see eq. 17). 
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Fig. 12. — Torque per unit disk mass on the planet as a function of normalized distance for the migrating and growing planets plotted in 
Figure 11. Left: Case with initial surface density at the initial orbit of the planet equal to E p 300gcm — 2 . Right: Case with initial surface 
density at the initial orbit of the planet equal to Ep fS 500gcm — 2 . The vertical scale is in units of GM a (M p /M s ) 2 /a, where a = a(t). The 
solid, long-dashed, dot-dashed, and short-dashed curves refer to times when M p = 6.0Me, 9.3 Me, 0.36 Mj, and l.OMj, respectively. 
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Fig. 13. — Migration with different initial conditions. Solid 
curve: Orbital radius evolution of a planet with initial mass 
M p = 5 that interacts with a three-dimensional disk hav- 
ing initial surface density at the planet's initial radial position 
E p Ri 300gcm~ 2 at ao = 5.2 AU (same as the upper migration 
track plotted in Fig. 11). It has mass M p = 0.3 Mj at a time of 
about 165 orbits (see Fig. 3, dashed line), when a ~ 0.92 ao. Dot- 
ted curve with solid circles: Orbital radius evolution of a planet 
with initial mass M p = 0.3 Mj that interacts with the same initial 
unperturbed disk density distribution as the solid curve case has 
at time t = 0. The planet starts at the same radius (a ~ 0.92 ao) 
as the solid curve where that planet has acquired a mass of 0.3 Afj. 
The difference in the two cases is that the solid curve case has 
a partially cleared gap when M p = 0.3 Mj (see Fig. 11, right 
panel), while the dotted curve case starts in a smooth unperturbed 
disk. Dashed curves: Orbital radius evolution of a planet accord- 
ing to Type I theory (eq. 17 and 18) for a planet of fixed mass 
M p = 0.3 Mj (lower curve of pair for saturated coorbital torques) 
and disk density at r = 0.92 ao for the unperturbed initial disk. 

The distribution is similar to the cases plotted in Fig- 
ure 1, although it is somewhat larger in magnitude, as 
expected by the lower sound speed of the gas and the 
two (rather than three) dimensions of the simulation. 
However, at later times the torque distribution changes 
character, with much larger values in the coorbital zone, 
within radial distances of about 2i?u — 0.1 a(t) from 
the orbit of the planet. In particular, there is substan- 
tial torque occurring in the radial band \r — a\ < Rn, 
where i?H = 0.046 a. We argued in section 2.2.2 that 
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Fig. 14. — Orbital migration of a Saturn-mass planet of fixed 
mass (M p = 0.3 Mj) in a cold (H/r = 0.03) and high mass disk 
(E p 670 g cm -2 at the planet's initial position). Mass is removed 
from the disk near the planet to prevent a mass buildup there. 
The planet is embedded in a two-dimensional disk and held on 
a fixed orbit for t r i a = 50 (dotted lines), 100 (solid line), and 
200 (long-dashed line) initial orbital periods. The dotted line with 
solid circles plots the migration track from a three-dimensional disk 
model with t r i s = 50 orbits. The orbital radius is in units of ao. 
For ao = 5.2 AU, the unit of time is 12 years. The predicted 
Type I migration tracks, assuming the planet does not open a gap, 
are plotted for a two-dimensional (short-dashed line) and a three- 
dimensional (short-dashed line with solid circles) disk. 

this region involves only coorbital torques (not Lindblad 
torques). We have verified that this torque is not origi- 
nating from within the planet's Hill sphere (see Fig. 27). 
The contribution from within the Hill sphere is about 
20% of the net torque at release time and generally less 
than about 10% at later times. 

Figure 15 shows that the planet is migrating on a 
shorter timescale than that of gap opening. At release, 
the planet is fairly symmetrically positioned in the gap. 
Later, the planet lies much closer to the inner edge of 
the gap than to the outer edge, and the gap is less deep. 
Such a situation would be expected to lead to slower or 
even outward migration according to Type I theory of 
Lindblad resonances, since the inner resonances (which 
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Fig. 15. — Axisymmetric radial density distribution E(r, t), of a 
disk containing a Saturn-mass planet (M p = 0.3 Mj), plotted as a 
function of radius at 3 times: the time of planet release f r i s = 100 
initial orbital periods (short- dashed line), t r i s + 10 initial orbital 
periods (long-dashed line), and t r i s + 20 initial orbital periods (solid 
line). The solid circles mark the planet orbital radii at these times, 
as the planet migrates inward. 
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Fig. 16. — Torque per unit disk mass on a Saturn-mass planet 
(M p = 0.3 Mj) in units of G M s (M p /M s ) 2 /a(t) as a function of 
the normalized distance from the planet, undergoing fast migra- 
tion, at 3 different times: the time of planet release t T \ s = 100 
initial orbital periods (short-dashed curve), t r \ s + 10 initial orbital 
periods (long-dashed curve), and t r i a + 20 initial orbital periods 
(solid curve). 



provide outward migration) are more strongly activated 
than the outer resonances (which provide inward migra- 
tion) due to the asymmetric density distribution near the 
planet. If fact, the slowing/stalling of inward migration 
due to the feedback from the inward disk density of a 
migrating planet was envisioned by Hourigan & Ward 
(1984) and Ward & Hourigan (1989) in their considera- 
tion of the inertial limit to planet migration. To quan- 
tify the effects of standard Type I torques, we apply the 
Type I torque distribution taken at the time of planet re- 
lease in Figure 16. In doing so, we are ignoring pressure 
effects on dT/dM due to the changing gap shape. We 
determine the Type I torques at the later times by in- 
tegrating this torque distribution (appropriately shifted 
to the instantaneous position of the planet) over the disk 
mass distributions in Figure 15. The torque is then given 



by 

Ti{t) = 2tt f^(x, t xls ) E(r, t) r dr (20) 

where x = (r—a)/a and a = a{t). We find that the result- 
ing Type I migration rates at times of 10 and 20 orbits 
after release are outward and equal to a = 2 x 10~ 4 cto £Iq 
and 4 x 10 _5 aorio: respectively Clearly, results from 
the simulation are not consistent with the expectations 
of the usual Type I migration theory. Instead, we claim 
the effects of the corotation resonances are critical for 
migration here. 

In the model by OL06, fast migration is due to torques 
caused by a density asymmetry in the coorbital region be- 
tween gas on the leading and trailing sides of the planet. 
The gas on the leading side of the planet is trapped and 
contains gas acquired at other radii, while the trailing 
side contains ambient gas near the planet. The con- 
trast between the trapped and ambient gas is limited by 
viscous diffusion. The trapped gas is in a quasi-steady 
advective-diffusive equilibrium. The density asymmetry 
and thus the torque is caused by the motion of the planet. 

To test this model, we analyzed streamlines in the coor- 
bital region in the frame comoving with the planet. We 
determine the streamlines in the simulations by follow- 
ing the motion of tracer particles that move with the 
velocity of the gas (see Appendix D.l). In Figure 17 
we plot coorbital streamlines near before the planet is 
released, i.e., while the planet was on a fixed station- 
ary orbit. The figure shows good agreement between 
the simulation and theory. The streamlines are symmet- 
ric between the leading (4> > 4>p) an d trailing (cp < 4> p ) 
sides of the planet. Figure 18 shows the streamlines af- 
ter the planet is released, while the planet is migrating. 
Strictly speaking, these are not streamlines in the sim- 
ulation case but trajectories, since the flow is not in a 
strict steady state in the comoving frame of the planet 
because the planet is migrating at a variable rate. The 
theoretical streamlines depend on the planet-to-star mass 
ratio and the migration rate of the planet. They are cal- 
culated assuming a steady state and constant migration 
rate by means of the linear perturbation model of OL06. 
The theoretical streamlines were calculated by using in- 
termediate parameter values from the simulation during 
the interval of planet migration: d = — 0.002 oo^o and 
r = 0.85 ao- The simulated and theoretical streamlines 
in Figure 18 are in approximate agreement. They show 
closed streamlines on the leading side of the planet's az- 
imuthal motion. They contain the trapped gas described 
above. The open streamlines on the trailing side of the 
planet involve ambient gas that streams outward past 
the planet. The smaller closed streamlines are centered 
at about the same azimuth in the two plots, about 0.2 it 
ahead of the planet. Figure 19 shows that the gas density 
asymmetry in the coorbital region between the leading 
and trailing sides of the planet increases with time. The 
unperturbed background density increases with time as 
the planet encounters higher density gas in its inward 
migration. Notice that the density increase is higher on 
the trailing side of the planet than on the leading side. 
This result suggests that the trapped gas approximately 
retains its initial density as the planet migrates. The gas 
on the trailing side more fully reflects the local density. 
The density asymmetry then gives rise to the dominant 
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Fig. 17. — Trajectories of gas in the coorbital region near a Saturn-mass planet (M p = 0.3 Mj) on a fixed circular orbit with orbital 
radius a = ao- The trajectories are determined in the comoving frame of the planet, which is located at the origin. Disks properties are 
the same as in Figure 14. The left panel shows the results of our simulation and the right panel shows results given by a theoretical model 
(OL06). 




0.0 




0.0 



Fig. 18. — Trajectories of gas in the coorbital region near a Saturn-mass planet undergoing fast migration (solid curve in Fig. 14). The 
trajectories are determined in the comoving frame of the planet, located at the origin. The thicker lines denote open trajectories that pass 
by the planet. The thinner lines denote closed trajectories containing trapped gas. Results from the simulation are presented in the left 
panel and results from a theoretical model (OL06) are displayed in the right panel. 



torque on the planet. 

The OL06 model does not determine the value of coor- 
bital corotational torque for a migrating planet. It does 
provide a detailed analysis for the noncoorbital corota- 
tional torque. In that case, the effect of migration is 
to amplify the standard coorbital torque for a nonmi- 
grating planet (Goldreich & Tremaine 1979). By anal- 
ogy, one might expect similar behavior in the coorbital 
case. The standard torque is proportional to the radial 
derivative of £(r)/_B(r), and hence of the disk vortensity 
—2B(r)/H(r), where S(r) is the disk axisymmetric sur- 
face density (azimuthally averaged surface density) and 
B{r) is the Oort constant. In the unperturbed disk model 
considered here, the vortensity is constant, and so the 
coorbital torque is predicted to be zero. However, the 
gap structure in the disk modifies the surface density 
S(r) near the planet (see Fig. 15), thereby providing a 
vortensity gradient. 

4.4. Conditions for Type III Migration 



The results in section 4.3.2 provide evidence for migra- 
tion dominated by coorbital torques, or Type III migra- 
tion. Within the framework of the OL06 model we gen- 
eralize the results of the simulations and describe some 
conditions that are favorable for Type III migration. 

As in section 4.3.2, we consider a planet of fixed mass. 
For a planet whose mass is large enough to open a gap, 
we apply the initial condition that the planet undergoes 
migration before steady gap formation completes, as was 
the case in section 4.3.2. 

We require that the planet does not strongly deplete 
the gas in the coorbital region as it migrates. This re- 
quirement implies that the migration timescale across the 
coorbital region be shorter than the timescale to clear a 
gap over that region. Figure 15 demonstrates that this 
condition holds for the model in section 4.3.2. 

To derive a crude estimate for this condition, we as- 
sume that the torques exerted by the planet lead to a 
local change in disk angular momentum over a region 
whose size is comparable to the coorbital region. Each 
one-sided torque (interior and exterior to the orbit of the 



16 



D'Angelo & Lubow 



2.5 






2.0 






1.5 






1.0 






0.5 


N / 




0.0 







-1.0 



-0.5 



(0. 



0.0 



0.5 



1.0 



Fig. 19. — Radial average of the surface density within the 
coorbital region (defined as having radial extent 2 Rh — 0.1 a(t) 
of the planet's orbit) as a function of azimuth at 3 different times: 
the time of planet release t T i s = 100 initial orbital periods (short- 
dashed line), t r i s + 10 initial orbital periods (long- dashed line), and 
t T i s + 20 initial orbital periods (solid line). 

planet) on the gas is capable of clearing a gap, while the 
net effect of both interior and exterior torques results in 
migration. The condition that the gap clearing timescale 
is longer than the migration timescale becomes 



M r > 



Mp 
A 



(21) 



where M c is the mass of the coorbital region and A is 
the dimcnsionlcss torque asymmetry 

\TJ - \TA 



A 



(22) 



where Tj and T e denote the torques interior and exterior 
to the planet's orbit, respectively. 

Another condition is that the migration rate be large 
enough that there a strong asymmetry in streamlines be- 
tween the leading and trailing sides of the planet, as 
seen in Figure 18. According to OL06, the asymme- 
try is strong for migration rates greater than \cla\ = 
1.45 57 a M p / M s . Since the migration begins as Type I 
migration (see short-dashed line in Fig. 16), the condi- 



tion is that 



«i 



> 



cla\- This condition is approximately 



V a 



(23) 



Notice that the condition is independent of planet mass, 
since both di and &a are linear in the planet mass. 

We now apply these conditions to the model simulated 
in section 4.3.2. The asymmetry parameter is estimated 
as A « 0.3, from the initial torque distribution in Fig- 
ure 16. Condition (21) is then satisfied for this model, 
since M c sw 7 M p . Condition (23) is also (marginally) 
satisfied, since the initial density is S p a 2 /M s = 2 x 10~ 3 
and (H/a) 2 = 9 x 10 -4 . None of the other models dis- 
cussed in this paper satisfy both conditions. 

It also appears that the condition that the planet mass 
is fixed (or slowly increasing) is important. The dashed 
curve on the right panel of Figure 27 in Appendix D sug- 
gests that a planet that grows in mass at the run-away 
rate would not undergo rapid migration long enough to 
move very far. The slow-down is partly due to gap open- 
ing that reduces the torques. 



The picture is then that a fixed mass planet, initially 
undergoing sufficiently fast Type I migration, develops 
strong coorbital torques due to asymmetric trapped gas. 
The situation is not simple, however. We saw in sec- 
tion 4.3.2, that the disk's feedback to the planet's mo- 
tion might slow or halt Type I migration, but the coor- 
bital torques allow the inward migration to continue. We 
found that the rate of the resulting migration is actually 
slower than Type I migration for a smooth disk. So the 
conditions (21) and (23) are suggestive only at this point 
and require further testing. 

5. SUMMARY AND DISCUSSION 

We have analyzed the evolution of migrating planets 
that undergo run-away gas accretion by means of multi- 
dimensional numerical simulations. The results agree 
with the predictions of Type I and Type II migration 
(see Figures 6 through 9) for a planet of time-varying 
mass that wc obtain from simulations. The set of sim- 
ulations include cases with disk densities as low as the 
minimum mass solar nebula value and as high as 5 times 
that value, viscosities a w 0.004 and 10 times that 
value (also 50 times that value, a f=a 0.2, in a test re- 
ported in Appendix B), disk temperatures corresponding 
to H/r = 0.05 and a colder case oiH/r = 0.04. The mass 
accretion rates onto the planet (see Figures 3 and 5) are 
in general agreement with previous determinations based 
on fixed mass planets on fixed circular orbits, when com- 
paring cases with the same planet mass and disk prop- 
erties for which the growth timescale is longer than the 
gap opening timescale. Planet mass growth rates can be 
understood in terms of accretion within the Bondi radius 
at lower planet masses, accretion within the Hill radius 
at intermediate masses, and accretion through the gap 
at higher planet mass (see Fig. 4). Mass growth rates 
typically peak at the intermediate planet masses, a few 
tenths of a Jupiter mass. 

An important diagnostic for the nature of the disk- 
planet torques is the scaled torque density distribution 
per unit disk mass as a function of the scaled radial dis- 
tance from the planet. In the linear regime of the stan- 
dard theory of disk-planet resonances, for a fixed form of 
the disk density distribution dlnE/dlnr and gas proper- 
ties (sound speed and viscosity), this scaled torque distri- 
bution should be universal, independent of planet mass 
and disk density value. We verified this universality for 
low mass planets, although the distribution varies some- 
what with planet mass for larger mass planets that open 
gaps (compare low mass cases in Figures 1,7, and 12). 

There is no fundamental distinction between the 
torques involved in Type I and Type II migration. This 
follows from the near independence of the scaled torque 
density distribution with planet mass. Previous concepts 
of Type II suggested that a planet in a clean gap would 
migrate inward like a test particle that follows the disk 
accretion. Here we describe a view for cases where the 
gap is not completely clear of material. The difference 
in Type I and Type II rates is due to the mass density 
distribution of the gas that multiplies torque density in 
determining the torque on the planet. In Type II mi- 
gration, the density distribution adjusts so that the net 
torque on the planet causes it to migrate at approxi- 
mately the viscous evolution rate of the disk. The transi- 
tion between the two forms of migration is quite smooth. 
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Fig. 20. — Orbital migration of a planet undergoing run-away 
gas accretion. Orbital radius in units of ao (5.2 AU), as a func- 
tion of time in units of the initial orbital period (rs 12 years). 
The initial planet mass is 5 Mg . The initial surface density is 
S p = 3 X 10 -4 M a ay 2 £S 100gcm~ 2 at the planet's initial or- 
bital radius and H/r = 0.05. Solid curve: Results from the same 
simulation as in Figure 6 for a migrating, mass-gaining planet. 
Short- dashed curves: Predictions based on Type I migration the- 
ory, obtained by solving equations (17) and (18), taking into ac- 
count the time variable mass of the planet (Fig. 3, solid line) and 
roughly accounting for the gas depletion near the planet (see text). 
The upper (lower) dashed curve is for migration with unsaturated 
(saturated) coorbital torques. 

To illustrate this point, we plot in Figure 20 the orbital 
evolution of the planet through the phase where gap for- 
mation sets in. The figure shows that the migration rate 
can be accounted for by standard Type I theory, cor- 
rected for the gas depletion in the gap region, although 
there is no unique prescription to do this. In the figure, 
the theoretical curves (dashed lines) are determined by 
Type I migration theory, equations (17) and (18), with 
the density S taken to be the average value in a radial 
band of half- width 0.15 a (a typical gap width) centered 
on the planet. 

For a given planet mass, the torque density diagnostic 
reveals that the distribution is not strongly affected by 
migration or accretion (Figure 1 is quite similar to Fig- 
ures 7 and 12 for similar planet masses). In particular, 
there is no evidence for strong coorbital torques. How- 
ever, in a certain case, we do find evidence for strong 
coorbital torques, or Type III migration. This case has a 
planet of fixed mass, 0.3A/j, that is immersed into a cold, 
smooth disk. The planet is held at a fixed orbit for rela- 
tively short time (~ 100 orbits) before being released, so 
that gap clearing is incomplete. The torque distribution 
at the time the planet is released follows the expecta- 
tions of the standard theory for nonmigrating planets. 
This can be seen by comparing the curve for the 0.3 Mj 
case in Figure 1 (dot-dashed curve) with the short-dashed 
curve (time t = t r i s ) in Figure 16. (There are differences 
in the magnitude of the distributions because of differ- 
ences in gas sound speeds and dimensionalities of the 
calculations, but the forms of the distributions are sim- 
ilar). At later times, Figure 16 reveals a transition to a 
completely different torque distribution, where coorbital 
torques play a critical role. 

The strong coorbital torque can be understood in terms 
of an asymmetry in the streamlines between the leading 
and trailing sides of the planet, in accord with the an- 



alytic model of OL06 (see Figures 17 and 18) and also 
along the lines of Artymowicz (2004). This asymmetry 
causes trapped material to persist on the leading side of 
the planet which has a different density from the ambi- 
ent gas that flows on the trailing side (see Fig. 19). The 
asymmetry gives rise to the coorbital torque. We suggest 
some criteria for this form of migration (see section 4.4). 
More exploration is needed to test them. 

Although we find evidence for Type III migration, the 
conditions required appear somewhat artificial, i.e., in- 
complete gap clearing (nonequilibrium gap) of a cool disk 
with a planet of fixed mass. It is not clear whether and/or 
how conditions for Type III could arise in a more plau- 
sible evolution scenario. 

We have generally assumed that the planet is able 
to accrete almost all gas the disk is able to provide 
(so-called run-away gas accretion). For a disk viscosity 
v > 1 x 10 -5 q,q flo, the mass accretion rates arc large. 
The time to build a 1 Mj planet starting with a 5 Me 
planet in a minimum mass solar nebula is shorter than 
~ 10 5 years, substantially less than the observationally 
determined disk lifetimes of ~ 10 6 years (Haisch et al. 
2001; Flaherty & Muzerolle 2008). Over this 10 5 year 
time interval, the planet has radially migrated inward by 
only ~ 20% of its initial radius. In other words, for these 
models, the mass doubling timescale for a M p < 1 Mj 
planet is short compared to the migration timescale and 
the disk lifetime. This situation stands in strong con- 
trast to the earlier phases of planet formation where the 
migration timescales are shorter than the planet mass 
doubling timescales and disk lifetimes (e.g., Ward 1997; 
Hubickyj et al. 2005). 

The run-away accretion rates pose some challenges 
for explaining the mass distribution of planets (Butler 
et al. 2006). Typical accretion rates in T Tauri stars 
are ~ 1 X 10 _8 M Q per year (Hartmann et al. 1998). 
For a steady-state unperturbed disk (without a planet), 
the accretion rate is given by ZitvYi. The initial disks 
considered in this paper are not in a steady state, but 
this accretion rate provides a reasonable estimate. For 
the minimum mass nebula model (Fig. 3), the kinematic 
turbulent viscosity we typically adopt, v ~ 10~ 5 G.gf2o, 
was chosen so that the accretion rate evaluates to about 
this same value, ~ 1 x 10~ 8 Mq per year. In the case of 
a planet embedded in a disk, if there is a comparable ac- 
cretion rate onto a planet of mass M p < lMj (as found 
by Lubow & DAngelo 2006), then the mass doubling 
timescale for a Jupiter-mass planet is about 10 5 years, 
consistent with what we found in the simulations in this 
paper. But then it is not at all clear why planets would 
not almost always achieve masses higher than ~ 1 Mj, 
in contradiction with the observed mass distribution of 
extra-solar planets and the case of Saturn. Special timing 
for disk dispersal could be invoked, but may be artificial. 

There are a few possible explanations. A colder disk 
(H/r < 0.05) would experience stronger tidal trunca- 
tion effects from a Jupiter-mass planet, as H becomes 
significantly smaller than i?H — 0.07 a. This effect could 
certainly reduce the accretion rate by a large enough fac- 
tor (say 10), so that the mass doubling timescale for a 
~ lMj planet would become of order 10 6 years, consis- 
tent with the suggestion by Dobbs-Dixon et al. (2007). 
The issue then is how cold the disk would need to be. 
The model in this paper with H/r = 0.04 (Fig. 5, long- 
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dashed line) has the same unperturbed overall disk ac- 
cretion rate as the H/r = 0.05 case (Figures 3 and 5, 
solid line), since v and the initial surface density S(r) 
are the same. The accretion rate onto the planet from 
the colder disk (H/r = 0.04) differs from the case with 
H/r = 0.05 by only 10%. Model h in Lubow & D'Angelo 
(2006) for a nonmigrating planet of fixed mass, having 
H/r = 0.03, has an accretion rate that is ~ 25% less 
than model b, which has H/r = 0.05 and the same un- 
perturbed overall disk accretion rate. Consequently, in 
this disk thickness range (0.03 < H/r < 0.05), we find 
that the reduction in the accretion rate onto the planet 
due to cooler disks is not significant. For higher mass 
planets (5-10 Mj), we expect that the accretion rate will 
be reduced by tidal truncation effects to a level where the 
planet mass doubling timescale is comparable to the disk 
lifetime, as found in previous studies of planets on fixed 
orbits (Lubow et al. 1999; Bate et al. 2003; D'Angelo 
et al. 2003). This effect may set the upper limit to planet 
masses. 

Another possibility is that there is a feedback effect 
that limits the gas accretion rate. Perhaps the heating 
of the protoplanet envelope by impacting solids contin- 
ues to later times than is assumed in the standard core 
accretion model. Depletion of disk solids near the planet 
occurs in the standard core accretion model, when planet 
migration is not included. With migration, it is possible 
that continued accretion of disk solids would occur (e.g., 
Alibert et al. 2005), resulting in continued heating that 
could limit the gas accretion rate further. It is not clear 
how well this possibility works, since planetesimals will 



get trapped into resonances as the planet migrates (Zhou 
& Lin 2007). Having a higher mass solid core is prob- 
lematic in the case of Jupiter, whose solid core mass is 
thought to be a small fraction of the total mass (see Guil- 
lot 2005 and discussion in Lissauer & Stevenson 2007). 
It is also possible that winds emanating from the cir- 
cumplanctary disk within the planet's Hill sphere could 
reduce the accretion rate onto the planet. Magnetically 
driven winds are believed to play an important role in 
the case of young stars (Blandford & Payne 1982; Pu- 
dritz & Norman 1986; Shu et al. 1994). There are likely 
differences in the flow properties from the stellar outflow 
case (Fcndt 2003). However, it is not clear that the winds 
would be able to expel a large enough fraction (say 90%) 
of the accreting gas to sufficiently reduce the accretion 
rate onto the planet. 
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APPENDIX 

A. NUMERICAL SENSITIVITY STUDY 

We conducted several tests to assess the sensitivity of the results presented in section 3 to the choice of various 
numerical parameters. Since the main objective of that section is the mass and orbital evolution of an embedded 
planet, we present here quantitative comparisons of planetary masses and orbital radii as function of time. 

A.l. A Resolution Test 

For purposes of a resolution study, we performed a three-dimensional calculation in which the grid resolution is 
raised by a factor of 3/2 over the standard resolution (see section 3) in each coordinate direction, throughout the 
entire disk domain, and on all grid levels. Note that such an increase implies an overall refinement gain of a factor 
(3/2) 3 ~ 3.4, in terms of volume resolution of the system or number of grid elements. Nested grids cover extended disk 
regions, so that the planet always remains in the domain described by the most refined grid during the calculation. We 
focus on the disk model with initial surface density at the planet's initial position £ p = 9 x 10 _4 M S <Zq 2 (300 gem -2 ), 
H/r = 0.05, and v = 1 x 10 -5 a^O . The planet's mass and orbital radius evolution, obtained at standard grid 
resolution, arc shown in Figure 3 (dashed curve) and Figure 11 (top-most solid curve), respectively. In order to carry 
out a quantitative comparison, results (for both M p and a) from the two calculations, which will be labelled as 1 and 
2, are averaged over time intervals of half of the (initial) orbital period and then relative differences are computed as 



AX 



Xi — x$ 



x I + x 2/ - < A1 > 

where X is either M p or a. In order to give time-averaged estimates of the relative differences, over the course the 
calculations, we perform a running-time average of quantity AX/X, which is defined by 

/ AA \ If* AX . , k . 

(ttW,-*' < A2 » 

Results are shown in the left panels of Figure 21. The largest differences are observed in the results for the mass 
evolution (top-left), after the onset of the rapid accretion phase (see dashed line in Fig. 3). The average difference, 
over the entire evolution, stays within 10-15%. The average relative difference between the evolution results of the 
orbital radii (bottom-left) is much smaller and remains well within 1%. 

In order to test whether the orbital evolution is consistent with Type I migration theory of Tanaka et al. (2002) at 
our standard resolution, we set up a three-dimensional disk model with a planet that grows at a prescribed rate and 
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Fig. 21. — Left: Running-time average, defined by equation (A2), of the relative differences between two three-dimensional calculations, 
whose numerical resolutions differ by a factor of 3.4 in terms of grid elements (see text for details). Most of the difference in the mass 
growth of the planet (top) is accumulated between about 100 and 150 orbits, during the early phases of the rapid mass growth, when 
M p ~ 1 X 10 -4 (see dashed curve in Fig. 3). Overall, the running-time average of AM P /M P is in the range 10-15%. Relative differences 
between the evolution of the orbital radii (bottom) remain negligible over the course of the calculations and the running-time average 
is contained within 1%. Right: Orbital migration of a planet that grows at a prescribed rate in a three-dimensional disk whose initial 
(unperturbed) surface density has slope s = — dlnT, p /dlna = 3/2. Results from a calculation (solid curve) are compared to predictions 
of Type I theory (dashed curve, see section 3.2.1). The top panel shows the orbital radius evolution, whereas the bottom panel shows the 
migration rate as a function of the planet mass. 

whose initial mass is M p = 1 x 10~ 5 M s (about 3 Me). The initial (unperturbed) surface density has slope s = 3/2, 
so that di(t) oc M p (t). In the right panels of Figure 21, results from the simulation (solid line) are compared against 
predictions of Type I theory (dashed line, see section 3.2.1) for both the orbital radius evolution (top panel) and the 
migration rate as a function of M p (bottom panel). 

A. 2. Boundary Condition Effects 

Boundary conditions may play some role and affect the late stages of the system's evolution, especially when M p 
becomes on the order of a Jupiter mass. In our situation, the major concerns that may arise are related to the 
positions of the grid radial boundaries. The finite extent of the inner radius of the grid (i? m in) may lead to an 
augmented depletion of the disk within the orbit of the planet. At the outer radial boundary (i? max ), reflection of 
waves may affect torques exerted on the planet. The first effect can be mitigated by reducing the inner radius of 
the grid or adopting nonreflective boundary conditions (e.g., Godon 1996, 1997), whereas the second can be largely 
suppressed by choosing i? max ^ a - However, while increasing the outer grid radius possibly lengthens the computing 
time only because a larger number of grid elements in the radial direction is required (for a given value of Ai?), 

3/2 

decreasing the inner grid radius directly affects the time-step of the calculation, which is proportional to R n [ in because 
of the stability criterion imposed by the Courant-Friedrichs-Lewy condition. 

The left panels of Figure 22 shows a comparison of results obtained from a three-dimensional calculation (T. p = 
9 x 10 -4 M s a^ 2 at t = 0, H/ r = 0.05, and v = 1 x 10" 5 a\ f2n) with nonreflective boundary conditions at i? m i n = 0.4 ao 
and one with i? m i n = 0.19 ao but outflow boundary conditions, as outlined in section 2.1.3. The position of the inner 
grid boundary affects the density distribution interior to the planet's orbit. The effect on the planet's mass is at the 
3% level, on average over the entire course of the simulation (top-left panel). The orbital radius evolution (bottom-left 
panel) displays average relative differences much smaller than 1%. In the right panels of Figure 22, results from a 
three-dimensional calculation with outer grid radius at i? m ax = 2.5 ao are compared to those from a calculation in 
which -Rmax = 4.9 ao- Both the evolution of planet mass (top-right panel) and orbital radius (bottom-right panel) are 
hardly affected by the position of the outer grid boundary. 

A. 3. Effects of Excluded Torgues 

The region of space in which material is gravitationally bound to the planet depends on several disk parameters 
(including H/r and v) and on the planet's mass. Calculations with fixed mass and fixed orbit planets indicate that, 
for H/r f=a 0.05 and v « 1 X 10~ 5 Oq Q,q, only material within about 0.3 i?H is gravitationally bound to the planet when 
5 Me < M p < 40 Me (Hubickyj et al. 2007). Analytical and numerical models of disk formation around a Jupiter- mass 
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Fig. 22. — Left: Running-time average of relative differences between results obtained from two three-dimensional models: one with inner 
grid boundary at radius i? m i n = 0.4 ao and nonrefiecting boundary conditions (Godon 1996, 1997) and the other with il m j n = 0.19 ao and 
outflow boundary conditions (see section 2.1.3). Average relative differences in planet's mass (top) and orbital radius evolution (bottom) 
are around 3% and <C~ 1%, respectively. Right: Same comparison between two three-dimensional calculations with outer grid boundary at 
radii fi max = 2.5 ao and 4.9 ao, respectively. Relative differences much smaller than 1% are observed in the evolution of both the planet's 
mass and orbital radius. 

planet suggest that such disks may extend over a distance of about Rr/A (or less) around the planet. In the calculations 
presented in section 3 and 4, torques originating within Ru/2 of the planet are not taken into account, which may 
include nonzero net torques exerted by unbound material. In order to estimate how this choice affects the evolution 
of the orbital radius, we also considered cases in which only torques from within 0.3 Rn are neglected. The models 
with initial surface density S p = 3 x 10~ 4 M s a^ 2 and 9 x 10~ 4 M s ag at r = a (H/r = 0.05, v = 1 X 10~ 5 a^flo), 
discussed in sections 3.1.2 and 3.2.2, are restarted at regular time intervals, covering entirely their respective mass 
range. The evolution is then integrated for time periods of > 100 orbits at each restart. We compare the evolution of 
orbital radii by measuring the relative differences Aa/a and find that for none of the cases considered |Aa/a| exceeds 
1%. We therefore conclude that excluded torques from unbound material have only marginal effects on the results 
presented in section 3.2 and 4.1. 



B. MIGRATION IN HIGH VISCOSITY DISKS 

Results presented in section 3.2 indicate that the orbital radius evolution of a growing planet can be described 
reasonably well in terms of standard Type I and Type II regimes of migration (as long as the local disk mass is 
comparable or larger than the planet's mass). This conclusion holds when the disk's kinematic viscosity is v ~ 
10~ 5 a 2 fio (see Fig. 6) as well as when v ~ 10~ 4 a 2 , f2o (see Fig. 9), which brackets a range of a-parameters between 
4 x 10~ 3 and 4 x 10~ 2 at the location of the planet. Here we present further analysis of a case with v = 1 x 10~ 4 a 2 , ilo 
and examine cases with v = 5x 10~ 4 a 2 , tto ~ 0.2), which still result in a form of Type I migration modified by the 
perturbed surface density of the disk. 

B.l. An Additional Model With v = 1 x 10~ 4 a 2 , fi 

As anticipated in section 3.2.2, some concern may arise when the viscous timescale, t v — r 2 /v, at i? max becomes 
comparable to the length of time over which the orbital evolution of the planet is calculated. However, it is unlikely 
that the disk's viscous evolution at R > i? max has a large impact on the results displayed in Figure 9 since t v at 
R = i?max is more than 6 times as long as the viscous timescale at a, about 10 4 (initial) orbital periods of the planet. 
Therefore, at this viscosity level, we may experience some effects only over simulations covering timescalcs longer than 
10 4 orbits (note that inward migration will increase even further this timescale). 

In order to address this issue more in detail, we set up a three-dimensional model with v = 1 x 10~ 4 agfirj and 
H/r = 0.05. We adopt the same parameters, numerical resolution (0.014 a of linear base resolution and 9 x 10~ 4 an of 
resolution in the coorbital region around the planet), and boundary conditions, as those of the simulations discussed 
in section 3. In this model, however, the grid radial boundary extends out to i? max = 6.7 ao- Therefore, the viscous 
timescale at i? max is a factor of at least 6.7 2 ~ 45 as long as t v at the orbital radius of the planet. Thus, we could 
in principle follow the planet's orbital evolution for tens of local viscous timescalcs. Additionally, to monitor the 
sensitivity of our results to boundary and initial conditions at i? m i n , we use the initial surface density represented 
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Fig. 23. — Left: Azimuthally averaged surface density of a three-dimensional disk with i/ = lx 10 -4 a§ f2o an d H/r = 0.05 at time t = 
(thinner solid line) and at times when M p = M p (t) is equal to 0.1 Mj (short- dashed line), 0.3 Mj (long-dashed line), and 1 Mj (thicker solid 
line). Right: Migration rate (da/dt) as a function of the planet mass (M p ) over the first 330 orbital periods of the simulation, during which 
time the planet's orbit is held fixed. The planet mass growth is prescribed. Migration rates are evaluated by means of Gauss perturbation 
equations. The upper and lower dashed curves indicate Type I migration rates predicted by equation (17) and (18), respectively. 




Fig. 24. — Left: Orbital radius evolution of a lMj planet in a disk whose kinematic viscosity is v = 1 X 10 _4 a§ f2 (H/r = 0.05). The 
solid line represents the model discussed in this Appendix while the dotted line with solid circles is the same as Figure 10. Note that the 
two models produce very similar migration tracks (see inset), although they use grids with different outer radial boundaries, respectively 
at 6.7 ao and 2.5 ao, and although the surface density profiles at r < a are different. In both calculations, the planet's orbit is held fixed 
over the first t r i s (initial) orbital periods (see text for details). The dashed (straight) line has a slope about equal to —7 X 10 — 5 aof2o. 
Right: Cumulative torque at time t = t T \ a + 600 initial orbital periods (thicker curve), in units of GM a M p /a, where a = a(t). See text for 
an explanation of the thinner curve. Nearly all torque is exerted by material within a radial distance of 0.25 a from the planet's orbit. 



as a thin solid line in the left panel of Figure 23. The initial mass density, p, is related to £ at t = as described 
in section 2.1.2. Material near the planet is removed from the disk, but its mass is not added to the planet's mass. 
Instead, for the purposes of this test, the planet mass is increased, at a prescribed rate, from M p = 1 x 10~ 5 M s (3 Me) 
to 1 x 10~ 3 M s (1 Mj) over about 330 orbital periods. The orbit of the planet is held fixed (a = ao) during this period 
of time. 

The left panel of Figure 23 shows the azimuthally averaged surface density at times when M p = 0.1 Mj (short- dashed 
line), 0.3 Mj (long-dashed line), and lMj (thick solid line). By means of Gauss perturbation equations (e.g., Beutler 
2005), we measure migration rates da/dt, which result from disk's gravitational forces (while a = ao), as a function 
of time and hence of M p = M p (t). The right panel of Figure 23 displays these static migration rates (solid line) 
compared to Type I rates (dashed curves) yielded by equation (17) and (18). The upper (lower) dashed curve refers 
to the unsaturated (saturated) coorbital corotation torques in the linear theory (see section 3.2.1). As observed in the 
figure (right panel), initial migration rates agree with those predicted by linear theory. The reduction of |d|, which 
peaks at M p /M s ~ 10~ 4 , is likely related to the onset of nonlinear effects (Masset et al. 2006). 

At time t r ] s = 340 orbital periods, the planet is released and allowed to change its orbit in response to disk's torques. 
We recall that, at this stage, the planet's mass is constant and equal to 1 Mj. The migration track is displayed as a 
solid line in Figure 24 (left panel). The planet's orbit is integrated for about 0.5 t v at ao- The dashed (straight) line 
has a slope about equal to —7 x 10~ 5 aof2o or —0.7v/a. For comparison purposes, also plotted in the left panel of 
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Fig. 25. — Left: Azimuthally averaged surface density of simulations with disk kinematic viscosity v = 5 X 10 -4 <2q Qo and H/r = 0.05. 
The planet mass grows from M p ~ 3 Me to lMj, at a prescribed rate, over roughly 330 initial orbital periods. The planet's orbit is fixed 
during this time interval. Thin solid line: Initial surface density distribution (same as the thin solid line in the left panel of Fig. 23). 
Thick solid line: Averaged surface density at time when M p = 1 Mj from a three-dimensional model whose outer radial boundary is 
fl max = 6.7 ao- Dotted line with solid circles: Averaged surface density when M p = lMj from a two-dimensional model with outer radial 
boundary at 13 ao. Dashed line: Case of a disk with M p = lMj and v = 1 X 10~ 4 aQf2o (same as the thick solid curve in the left 
•panel of Fig. 23). Right: Orbital radius evolution after release, t T \ s = 340 initial orbits, obtained from the three-dimensional model with 
Kmai = 6.7ao (solid line) and the two-dimensional model (dotted line with solid circles). Note that the orbital evolution covers more than 
1.8 local viscous timescales. The inset shows migration tracks obtained from the three-dimensional calculations over t„/2 viscous timescales 
at the initial orbital radius of the planet. The solid diamonds represent data from the case with i? max = 13 ao- The solid line represents 
the same model as in the main panel. 

Figure 24, as a dotted line with solid circles, are results shown in the right panel of Figure 9 (dotted line with solid 
circles) and obtained with the model discussed towards the end of section 3.2.2, which has a different disk density 
profile inside the planet's orbit and a different radial coverage of the disk. Despite these differences, the two models 
produce consistent and very similar migration tracks, indicating that disk regions at radii much smaller and much 
larger than the planet's orbital radius are not playing a determinant role. 

The thick curve in the right panel of Figure 24 plots the cumulative torque (in units of GM s M p /a), i.e., the torque 
per unit radius dT/dr integrated outward over radius, at t = t T \ s + 600 orbits. As also noted for the lMj case in 
Figure 2, almost all the torque is due to material within a radial distance of 0.25 a from the orbit of the planet. 
Figure 24 (left panel) indicates that the migration rate decreases over the course of the simulation. We find that this 
is not caused by a changing character of the torque per unit disk mass, dT/dM, but rather by a changing (azimuthally 
averaged) surface density profile around r ~ a. We calculate dT/dM at release time and the averaged surface density 
600 orbits after release. We then use equation (20) to obtain the expected cumulative torque, at time t = t T \ s + 600 
orbits, under the assumption that the intrinsic character of dT/dM remains unchanged over time. This is plotted as 
a thin curve in the right panel of Figure 24, along with the actual cumulative torque (thick curve). The difference 
between the two curves is less than f0%. 

B.2. A Model With i/ = 5x f 0~ 4 a§ fl 

We wish to examine here whether the migration trend observed in raising the viscosity from v ~ 10 -5 a\ f2 to 
v ~ f 0~ 4 <2q Oo persists at larger viscosity. As shown in the left panel of Figure 23 (solid line), when v = 1 x 10 -4 a\ f2o 
a Jupiter-mass planet is able to open only a shallow gap along its orbit (there is a drop in density of about a factor 3 
relative to the value just outside the gap). This is because gap opening conditions are not satisfied (see section 3.2.2). 
At larger disk viscosity we therefore expect an even shallower gap. 

We perform two three-dimensional simulations with the same setup as that outlined above but with v = 5 x 10 -4 a§ f?o 
(a = 0.2 at r = a ) and outer radial boundaries at R maK = 6.7 a and 13 ao, respectively. The linear base resolution 
is AR = a Ad = a A<fi = 0.014 a while the resolution in the coorbital region around the planet is 9 x 10~ 4 ao- The 
viscous diffusion timescalc, t v , at r = ao is approximately 320 (initial) orbital periods whereas t v at -R max = 6.7 ao is 
over 1.4 x 10 4 orbits (and 5.4 x 10 4 orbits at R max = 13 ao). We also consider a two-dimensional version of such models, 
having the same grid structure and resolution in the r-<p plane, and outer grid boundary located at r max = 13 ao (nearly 
68 AU from the central star). The planet mass grows, at a prescribed rate, from M p ~ 3 Me to 1 Mj over about 330 
periods (which is similar to t v at the planet position), while the planet's orbit is held fixed. The left panel of Figure 25 
shows the initial surface density (thin solid line) and the azimuthally averaged density profile when M p = 1 Mj for the 
three-dimensional model with i? max = 6.7 ao (thick solid line) and the two-dimensional model (dotted line with solid 
circles). As reference, the azimuthally averaged density for the case with viscosity ^=lxl0~ 4 aQf2o and M p = 1 Mj 
is also plotted as a dashed line (same as the thick solid curve in the left panel of Fig. 24). In the overlapping disk 
region, two- and three-dimensional calculations give consistent results. No significant deviations from the initial density 
distribution are observed at r > a. 

The viscosity condition for gap opening requires that M p /M s > 0.02 (see section 3.2.2), or M p > 20 A/j, in order for 
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Fig. 26. — Orbital evolution under the same disk conditions as for models in Figure 11. But, the planet's angular speed is imposed 
to be equal to the instantaneous keplerian value while the planet migrates in response to the nonaxisymmetric disk forces. Solid curves: 
Simulation results for orbital migration of a planet in a disk with initial surface density equal to S p = 9 X 10 -4 M s d^ 2 , or about 
300gcm~ 2 at ao = 5.2 AU (upper migration track), and E p = 1.5 X 10 -3 M a a^" 2 , or about 500gcm~ 2 (lower migration track). Dashed 
curves: Predicted orbital migration according to Type I theory, equations (17) {upper curve of pair for unsaturated coorbital torques) and 
(18) (lower curve of pair for saturated coorbital torques). Analogous calculations executed for density distributions and disk thicknesses 
used in section 3 produce migration tracks that differ by ~ 1%, over the entire planet mass range, from those displayed in the left panels 
of Figures 6, 8, and 9. 

gravitational torques to overcome viscous torques. In fact, the density distribution in the left panel of Figure 25 [thick 
solid line and dotted line with solid circles) shows a form of rather shallow gap. Hence, Type II migration should not 
be expected. 

The planet is released from its fixed orbit at time t r \ s — 340 (M p — 1 Mj for t > t T i s ) and the orbit is integrated for 
1.8 t v viscous timescales at r = ao- Figure 25, [right panel) displays the orbital radius evolution for the two-dimensional 
[dotted line with solid circles) and three-dimensional [thick solid line) simulations. The two migration tracks in the 
main panel closely follow one other. A comparison between the results obtained from the three-dimensional models, 
over i„/2 at the initial orbital radius of the planet, is shown in the inset. Again, there is no indication that disk's 
evolution at r > o has a significant influence on planet's migration. As for the case discussed above, nearly all the 
torque is accumulated by material within a radial band |7" — a.| < 0.25 a centered on the planet's orbit. 

The rate of migration after release time is approximately —8 x 10~ 5 ao Oo, which is similar to that of the solid curve 
in the left panel of Figure 24. This near equality is expected for Type I migration, since it is independent of the level 
of disk viscosity. We use the torque per unit disk mass, dT/dM, at t = t r \ s from the model with v = 1 x 10 _4 aQ Hq 
discussed above in Appendix B.l and the averaged surface density profile in the left panel of Figure 25 [thick solid 
line). By applying equation (20), we estimate the total torque expected under the assumption that dT/dM has 
similar shapes in the two models. There will be some dependence of dT[r)/dM on viscosity, since viscosity affects the 
resonance widths. This estimate yields a migration rate that agrees within a factor of 1.7 with the value stated above, 
indicating that the intrinsic character of the torque per unit disk mass is roughly similar in these two cases. 

C. CORRECTIONS FOR DISK GRAVITY 

As discussed in section 4.1, there is a possible artificial torque that can act on a planet surrounded by a massive disk, 
when the planet responds to the gravity of the disk but the disk self-gravity is not included (Pierens & Hure 2005). 
This torque is a consequence of the disk's axisymmctric gravitational force in changing the planet's orbital rotation 
rate, but not changing the disk's rotation rate (since the disk is not self-gravitating). This artificial difference leads to 
a shift in disk resonances that in turn leads to an artificial increase in the planet's inward migration. It can largely be 
remedied by forcing the planet to rotate at the local keplerian rate, i.e., at the same speed as the gas rotates apart from 
effects of gas pressure. In this prescription, the planet responds to the nonaxisymmetric forces of the disk that result in 
migration, while undergoing orbital motion at the keplerian rate. This scheme is in reasonable accord with simulations 
that include the full effects of disk self-gravity (Baruteau & Masset 2008). The full effects of self-gravity cause a slightly 
faster migration rate than this approximation suggests. We have carried out three-dimensional simulations with such 
imposed keplerian planetary orbits for various disk mass cases discussed in sections 3 and 4. For mass distributions 
and disk thicknesses, as those applied in section 3, migration tracks show negligible differences, over the entire planet 
mass range. The only cases that produce changes beyond a few percent in migration rates are those in Section 4. In 
Figure 26 we plot the resulting migration for the same disk models as in Figure 11. The migration rates are slower, as 
found by Baruteau & Masset (2008). But, they arc still in approximate agreement with the predictions of migration 



24 



D'Angelo & Lubow 




Fig. 27. — Orbital migration of a 0.3 Mj planet in a cold and massive disk (see section 4.3 for details). The planet's orbit is held fixed 
for i r i s initial orbital periods. Left: Comparison between two-dimensional models with (solid lines) and without (dotted lines) accreting 
boundary conditions near the fixed mass planet (t T \ a = 100). The curves marked with solid circles represent migration tracks obtained 
by excluding torques from within the planet's Hill sphere. Right: Migration tracks from three-dimensional models with sudden release 
(*rls = l)i an d fixed (solid and dotted lines) and variable (dashed line) mass planets. The solid and dotted curves are for an accreting and 
nonaccrcting planet, respectively. The long-dashed curve represents a case in which the initial planet mass (0.3 Mj) is augmented by the 
mass of the gas within -Rh/4 of the planet. 



theory. 

D. ADDITIONAL TESTS ON FAST MIGRATION 

Simulations of the orbital evolution of a fixed Saturn-mass planet [M p = 0.3 Mj) in a cold [H/r = 0.03) and massive 
disk (Sp = 2 x 10~ 3 M s Oq 2 w 670 gem -2 at the planet's initial orbital radius) can lead to a buildup of gas within the 
planet's Hill sphere, which is eventually halted when a sufficiently large pressure gradient is established. The mass of 
material that accumulates around the planet can exceed the planet's mass, with possible effects on migration rates. 
In order to prevent the accumulation of gas within the Hill sphere, in the models presented in section 4.3, we applied 
accreting boundary conditions near the planet, without adding the gas mass to the planet mass. In this Appendix 
we wish to reconsider the nonaccrcting configuration (as in MP03 and DBL05). The nonaccrcting approach may be 
considered to be crudely simulating a case where some process prevents the planet from gaining further mass. 

In the left panel of Figure 27, migration tracks from a two-dimensional model with an accreting planet [solid curves) 
are compared to those obtained from a two-dimensional model with a nonaccrcting planet [dotted curves). The gas 
masses within the Hill spheres are drastically different: ~ 0.07 M p and ~ 1.6 M p in the accreting and nonaccrcting 
planet cases, respectively. In the nonaccreting case, the mass of the gas (~ 1 M p ) within the bound region (see Fig. 29, 
right panel) should be added to the inertial mass of the planet, thereby slowing migration (Papaloizou ct al. 2007). 
This effect is indeed seen at early times (less than 10 orbits after release). The migration of the nonaccreting planet 
with bound gas [dotted curve) is slowed by about a factor of 2 relative to the accreting case [solid curve). At later 
times, the migration rates are closer, although it is not clear why. The reason may be related to our determination 
that the mass of "bound" gas decreases at later times in the nonaccreting case (see also right panel of Fig. 29). 

The left panel of Figure 27 also plots the orbital evolution when torques from within the Hill sphere arc not taken 
into account [solid and dotted curves with solid circles). The similarity of these migration tracks to the other plotted 
in the figure indicates that torques from gas in this region do not dominate the migration rates in this particular case. 

In the nonaccreting case, dense gas that accumulates around the planet could be thought of as forming an envelope, 
once it becomes bound to the planet. A massive envelope would then participate in both the gravitational and 
inertial mass of the planet. We set up a three-dimensional model, with the same disk properties mentioned above (see 
section 4.3.1 for details), and planet mass M p = M c + M e , where M c = 0.3 Mj is a "core" mass and M e = M e [t) 
is the mass of the gas within i?H/4 of the planet. Given the large initial mass in the coorbital region (~ 2M,j), the 
planet rapidly gains mass, growing beyond 1 Mj in less than 25 initial orbits. The planet is released in a smooth disk 
after a few orbits. The orbital radius evolution is shown as a dashed line in the right panel of Figure 27, together with 
those from three-dimensional models with a fixed mass planet and accreting [solid line) and nonaccreting [dotted line) 
boundary conditions near the planet. The initial migration rates are similar in all three configurations but migration 
starts to rapidly slow down when the planet mass, in dashcd-line case, grows beyond M p ss 0.8 Mj. This behavior 
resembles that seen in Figure 13 [dotted line with solid circles). 

Figure 28 displays numerical convergence tests for the accreting [left) and nonaccreting [right) planet models pre- 
sented in the left panel of Figure 27. The two simulations in the left panel have coarsest (linear) resolutions in the 
coorbital region that differ by a factor of 4, in both radial and azimuthal directions. Calculations in the right panel 
have resolutions in the coorbital region around the planet that differ by a factor of 2 in each direction. 
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Fig. 28. — Numerical convergence tests for two-dimensional models with and without accreting boundary conditions near the planet (see 
Fig. 27). The planet's orbit remains fixed for f r i s = 100 initial orbital periods. Left: Comparison between simulations with an accreting 
planet and coarsest grid resolution in the coorbital region Ar = ao A</> = 0.014 ao (dashed curve ) and 3.5 x 10~ 3 (solid curve). Right: 
Comparison between models with a nonaccrcting planet and a linear resolution in the coorbital region around the planet of 0.02 Rn (dotted 
curve) and 0.01 Ru (dashed curve). 
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Fig. 29. — Tracer particles deployed within about 2 i?n of a 0.3 Mj planet orbiting in the cold disk model (H/r = 0.03) discussed in 
this Appendix and in section 4.3. The plot shows the distance, in units of ijjj, from an accreting (left) and nonaccreting (right) planet. 
Distances of trajectories that return to the disk are indicated as thin lines, otherwise they are indicated as a thick lines. 

D.l. Gas Bound to the Planet 

In the calculations with an accreting planet, torque contributions from within Rn/2 of the planet are ignored. By 
following fluid paths, here we show that most of this material is captured and eventually accreted by the planet. 

In a nonstationary flow, streamlines can be used as a proxy for fluid trajectories only over short distances and periods 
of time. Therefore, we track trajectories of fluid parcels by deploying tracer (massless) particles in the flow and then 
following their motion. This procedure allows us to obtain a reliable determination of fluid paths regardless of whether 
the flow is close or far from steady state. 

The equations of motion of each particle are integrated every hydro-dynamical time-step by interpolating the velocity 
field at the particle's location and by advancing its position in time via a second-order Rungc-Kutta method. Both 
spatial and temporal interpolations are performed by using the velocity field with the highest resolution available, i.e., 
that belonging to the most refined grid level in which the particle resides. The spatial interpolation is based on a 
monotonized harmonic mean (van Leer 1977), which is second-order accurate and capable of handling discontinuities 
and shock conditions. Hence, trajectories are formally second-order accurate in both space and time. 

Here we employ tracer particles to estimate the size to the region occupied by gas bound to nonmigrating planets. 
Tracers are deployed in the disk within about 2i?H of a Saturn-mass planet (M p =0.3 Mj). We use both models with 
and without accreting boundary conditions near the planet discussed in this Appendix (Fig. 27, left panel) and in 
section 4.3. Figure 29 shows the distance from the planet, S, of particles as a function of time. Tracers are deployed at 
a time t^, when the mass within the Hill sphere has reached a nearly steady value. The distance along the trajectories 
is normalized to the Hill radius, i?H- The left panel refers to the accreting planet case, whereas the right panel refers 
to the nonaccreting planet case. Thin curves mark distances of trajectories that are not captured and thus return 
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to the disk. Thick curves mark distances of trajectories that are captured in the planet's gravitational potential (left 
panel) or otherwise remain within Rn/2 of the planet, over the simulated evolution (right panel). In the accreting case, 
bound trajectories rapidly decay towards the planet and it is therefore possible to make a clear distinction between 
bound and unbound trajectories. In the nonaccrcting case, the distinction is less clear and may apply only over a 
given amount of time. 

REFERENCES 



Alibcrt, Y., Mordasini, C, Bcnz, W., & Winisdocrffcr, C. 2005, 

A&A, 434, 343 18 
Artymowicz, P. 1993, ApJ, 419, 155 3 

Artymowicz, P. 2004, in KITP Conference: Planet Formation: 

Terrestrial and Extra Solar 8, 17 
Baruteau, C, & Masset, F. 2008, ApJ, 678, 483 10, 23 
Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, 

MNRAS, 341, 213 1, 2, 3, 7, 18 
Bcutlcr, G. 2005, Methods of celestial mechanics. Vol. I: Physical, 
mathematical, and numerical principles (Methods of celestial 
mechanics. Vol. I / Gerhard Beutler. In cooperation with Leos 
Mervart and Andreas Verdun. Astronomy and Astrophysics 
Library. Berlin: Springer, ISBN 3-540-40749-9, 2005, XVI, 464 
pp. 99 figures, 11 in color, 32 tables and a CD-ROM.) 21 
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 18 
Bodcnhcimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 1 
Brydcn, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, 

J. C. B. 1999, ApJ, 514, 344 7 
Butler, R. P., Wright, J. T., Marcy, G. W., Fischer, D. A., Vogt, 
S. S., Tinney, C. G., Jones, H. R. A., Carter, B. D., Johnson, 
J. A., McCarthy, C, & Penny, A. J. 2006, ApJ, 646, 505 2, 17 
D'Angelo, C, Bate, M. R., & Lubow, S. H. 2005, MNRAS, 358, 
316 (DBL05) 1 

D'Angelo, G., Henning, T., & Kley, W. 2002, A&A, 385, 647 2 
D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 1. 2, 
7, 18 

D'Angelo, G., Lubow, S. H., & Bate, M. R. 2006, ApJ, 652, 1698 
1 

Dobbs-Dixon, I., Li, S. L., & Lin, D. N. C. 2007, ApJ, 660, 791 2, 
17 

Egglcton, P. P. 1983, ApJ, 268, 368 7 
Fendt, C. 2003, A&A, 411, 623 18 

Flaherty, K. M., & Muzerollc, J. 2008, AJ, 135, 966 17 
Godon, P. 1996, MNRAS, 282, 1107 19, 20 
— . 1997, ApJ, 480, 329 19, 20 

Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 15 
— . 1980, ApJ, 241, 425 2, 3 

Guillot, T. 2005, Annual Review of Earth and Planetary Sciences, 
33 493 18 

Haisch, K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 17 
Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 
495, 385 17 

Hourigan, K., & Ward, W. R. 1984, Icarus, 60, 29 14 
Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 
415 1, 5, 17 

Hubickyj, 0., Lissauer, J. J., D'Angelo, G., & Bodcnhcimer, P. 

2007, AGU Fall Meeting Abstracts, A6 19 
Kley, W. 1999, MNRAS, 303, 696 7 

Li, H., Li, S., Koller, J., Wendroff, B. B., Liska, R., Orban, C. M., 
Liang, E. P. T., & Lin, D. N. C. 2005, ApJ, 624, 1003 1 



Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 1, 8 

Lin, D. N. C, & Papaloizou, J. C. B. 1993, in Protostars and 

Planets III, 749-835 2, 9 
Lissauer, J. J., & Stevenson, D. J. 2007, in Protostars and Planets 

V, ed. B. Reipurth, D. Jewitt, & K. Keil, 591-606 18 
Lubow, S. H., & D'Angelo, G. 2006, ApJ, 641, 526 7, 17, 18 
Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 

2, 3, 18 

Marcy, G., Butler, R. P., Fischer, D., Vogt, S., Wright, J. T., 
Tinncy, C. C, & Jones, H. R. A. 2005, Progress of Theoretical 
Physics Supplement, 158, 24 2 

Masset, F. S., D'Angelo, G., & Kley, W. 2006, ApJ, 652, 730 8, 
21 

Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494 (MP03) 
1 

Mihalas, D., & Weibel Mihalas, B. 1999, Foundations of radiation 

hydrodynamics (New York: Dover, 1999) 2 
Nelson, A. F., & Benz, W. 2003, ApJ, 589, 578 1 
Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, 

MNRAS, 318, 18 1 
Ogilvie, G. I., & Lubow, S. H. 2006, MNRAS, 370, 784 (OL06) 1 
Paczyhski, B. 1971, ARA&A, 9, 183 7 

Papaloizou, J. C. B., Nelson, R. P., Kley, W., Masset, F. S., 
& Artymowicz, P. 2007, in Protostars and Planets V, ed. 
B. Reipurth, D. Jewitt, & K. Keil, 655-668 12, 24 

Pierens, A., & Hure, J.-M. 2005, A&A, 433, L37 10, 23 

Pollack, J. B., Hubickyj, O., Bodcnhcimer, P., Lissauer, J. J., 
Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62 1, 5 

Pudritz, R. E., & Norman, C. A. 1986, ApJ, 301, 571 18 

Shu, F., Najita, J., Ostriker, E., Wilkin, F., Ruden, S., & Lizano, 

S. 1994, ApJ, 429, 781 18 
Tanaka, H., Takeuchi, T., & Ward, W. 2002, ApJ, 565, 1257 4, 8, 

18 

Tanigawa, T., & Watanabc, S. 2002, ApJ, 580, 506 6 
van Leer, B. 1977, JCP, 23, 276 25 
Ward, W. 1997, Icarus, 126, 261 2, 8, 17 
Ward, W. R. 1986, Icarus, 67, 164 3, 4 

Ward, W. R. 1992, in Lunar and Planetary Institute Conference 
Abstracts, Vol. 23, Lunar and Planetary Institute Conference 
Abstracts, 1491-1492 4 

Ward, W. R., & Hourigan, K. 1989, ApJ, 347, 490 14 

Wuchterl, G. 1991, Icarus, 91, 53 1 

— . 1993, Icarus, 106, 323 5 

Yuan, C, & Cassen, P. 1994, ApJ, 437, 338 4 

Zhou, J.-L., & Lin, D. N. C. 2007, ApJ, 666, 447 18 

Zieglcr, U., & Yorkc, H. W. 1997, Computer Physics 
Communications, 101, 54 2 



